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PREFACE 


From  1971  through  1973,  a  new  sampled-data  processing  technl'  2  tor 
digital  signals  subject  to  colored  multiplicative  noise  was  developed 
and  subsequently  patented  by  the  Principal  Investigator,  at  NASA  Langley 
Research  Center.  In  1974,  a  contract  was  issued  by  the  Air  Force  Avionics 
Laboratory  to  determine  if  the  same  technique  which  provided  processing 
gain  against  diffuse  Doppler-spread  multipath  perturbations  could  be 
applied  to  anti -jam  processing. 

Anti-jam  processing  algorithms  were  produced  under  the  1974  contract, 
as  well  as  a  Monte  Carlo  simulation  package  for  performance  evaluation. 
Between  1976  and  1978,  substantial  evaluation  of  the  algorithms  was  per¬ 
formed  and  documented,  under  an  extension  of  the  contract. 

A  final  extension  of  the  contract,  through  April  1979  served  to  sup¬ 
port  investigation  of  means  for  implementing  carrier  phase  estimation  and 
bit  synchronization  with  the  detection  algorithms.  This  report  documents 
those  results  and  gives  recommendations^ for  further  research. 
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SECTION  I 


INTRODUCTION 

This  report  documents  further  research  under  the  subject  contract 
whose  previous  results  have  been  reported  in  [1,4].  The  basic  te  .inical 
problem  is  that  of  optimum  discrete-time  recursive  detection  of  binary 
signals  subject  to  additive  colored  and  white  noise.  Previous  results 
showed  that  the  Minimum  Probability  of  Error  detector  is  one  which  tracks 
the  colored  noise  and  subtracts  it  from  the  received  data.  The  related 
question  of  identification  of  the  statistics  of  the  colored  interfering 
process  was  extensively  investigated  in  Reference  1. 

The  research  effort,  documented  herein,  was  pointed  toward  several 
related  questions.  First,  it  was  desired  to  investigate  the  problem  of 
simultaneous  estimation  of  the  carrier  phase  references  required  by  the 
coherent  detection  algorithm.  It  was  desired  to  specifically  determine 
the  method  for  measuring  phase  and  also  the  augmentation  of  the  detection 
algorithm  required  to  operate  with  imperfect  phase  estimates. 

Next,  it  was  desired  to  investigate  the  possibility  of  non-coherent 
detection  with  interference  tracking,  with  application  to  Frequency-Shift- 
Keying  and  Differential  Phase-Shift-Keying. 

A  third  area  of  interest  was  to  determine  a  method  for  obtaining  bit 
synchronization  for  the  interference-tracking  detection  algorithms.  This 
would  then  lead  to  assembly  of  a  complete  algorithm  for  the  so-called 
IDEI  (Integrated  Detection,  Estimation,  Identification)  receiver. 

Finally,  it  was  desired  to  obtain  Monte  Carlo  evaluation  of  the  aug¬ 
mented  detector,  operating  in  an  environment  of  colored  plus  white  addi¬ 
tive  noise. 

All  of  the  desired  areas  are  investigated  below.  An  expected  result 
is  that  the  coherent  detector  performance  is  degraded  when  carrier  phase 
is  estimated  from  the  received  data.  An  unexpected  result  is  that  a  non¬ 
coherent  version  of  the  interference-tracking  detection  algorithm  does 
not  exist. 

Recommendations  are  given  on  further  research  which  may  lead  to  im¬ 
proved  performance  of  the  complete  IDEI  receiver. 
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SECTION  II 

COHERENT  DETECTION  WITH  PHASE  ESTIMATION 


1.  SIGNAL  AND  CHANNEL  MODEL 

Figure  1  shows  the  overall  model  of  the  signal  channel  and  signal 
processor.  A  conti nuous-time  signal,  a(t,m),  is  transmitted  through  the 
channel . 

4(t,m)  =  A(t;m)cos[u>ct  +  <f>(t;m)]  (1) 

In  (1),  A(  )  and  <j>(  )  are  the  envelope  and  phase  functions,  respectively, 
m  denotes  a  digital  symbol,  which  in  the  present  work  is  restricted  to 
the  binary  alphabet,  {0,1}.  Any  arbitrary  signal  waveform  may  be  repre¬ 
sented  in  the  form  of  (1). 

The  signal  is  subjected  to  additive  colored  and  white  noise,  as  per 
the  figure.  Then  the  bandpass  signal  plus  noise  process  is  translated  to 
baseband  in  two  separate  channels,  using  coherent  product  detection  with 

sinusoidal  reference  signals  which  are  in  phase  and  in  phase  quadrature  with 
the  unmodulated  carrier  signal.  Following  the  I-Q  demodulation,  the  two 
low-pass  signal  components  of  the  signal  vector  are  sampled  to  produce 
a  discrete-time  vector.  The  discrete-time  signal  is  then  processed  further 

/V 

to  recover  the  message  symbol  decisions,  m, 

The  I-Q  product  demodulators  require  reference  sinusoids  having  pre¬ 
cise  phase  references,  matched  to  the  phase  (zero)  of  the  unmodulated 
carrier  signal.  Since  this  phase  is  A  Priori  unknown,  the  phase  reference 
must  be  provided  by  the  signal  processor,  itself,  by  phase  estimation 
from  the  received  data  vector.  The  reference  phase,  so  produced,  is 
generally  a  function  of  time,  q>Q(t) ,  as  shown  in  Figure  2. 

Since  the  signal  phase  is  A  Priori  unknown,  the  signal  model  of  (1) 

may  be  augmented  with  a  random  (or  stochastic)  phase  term  <p  as 

•6 

a(t,m)  =  A(t;m)cos[oo  t  +  <j>(t;m)  +  $  ] 

C  4 

=  ^ (t;m)cosu)ct  -  4q(t;m)sinuct  (2) 

where 

a.(t;m)  =  A(t;m)[cos4>(t;m)cos<t>  -  sin<j>(t;m)sin4>  ] 

I  4  4 

<s_(t;m)  =  A(t;m)[sin<j>(t;m)cos<J>.  +  cos4>(t;m)sin4>.  J  (3) 

H  4  4 
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Figure  1.  Physical  Channel  and  Receiver  Models 


LPF 


s(t;m)  ♦  n(t) 


2  c<*[<uct  +  <£0(t)] 


-2  sin[cuct  +  ^(t)] 


Figure  2.  I-Q  Carrier  Demodulator 

are  the  in-phase  and  quadrature  low-pass  components  of  the  band-pass 
a (t;m). 

The  I-Q  components  of  4(t;m)  form  a  vector 


4.j(t;m) 

*q(t;m) 


cos<{)  -si n<p  I  A(t;m)cos<|>(t;m) 

*  A  =  a(t;m) 

sin<b  cose}),  A(t;m)sin<j>(t;m) 

■6  4  L-  _ 


Likewise,  the  additive  colored  and  white  noises  may  be  written  in  terms 
of  I-Q  components  as 


y,(t) 

*(t)  =  1 

yq(t) 


n.(t) 

;  n(t)  = 

"q(t> 


where  jf(t)  is  the  low-pass  I-Q  colored  interference  vector  and  ji(k)  is 
the  I-Q  data  vector,  z{t)  may  then  be  written  as 

z{t)  =  a(t;m)  +  y(t)  +  n(t)  (6) 


The  problem  of  detecting  the  digital  symbol,  m,  in  the  presence 
of  colored  noise,  white  noise,  and  unknown  signal  phase  is  essentially 
the  problem  of  processing  _z ( t )  to  make  an  optimum  decision  on  m.  inis 
problem  is  analyzed  in  some  detail  below. 


2.  JOINT  DETECTION  WITH  PHASE  ESTIMATION 

It  is  desired  to  reformulate  the  discrete-time  recursive  detection 
problem  of  [1]  for  the  present  case  where  the  signal  phase  is  unknown  and 
time-varying.  At  this  point  it  is  still  assumed  that  the  symbol  epoch, 
or  timing,  is  known.  The  decision  problem  is  based  on  processing  the 
discretized  I-Q  data  vector  of  (6).  That  is,  a  sequence  of  samples,  z(tk) 
is  processed  recursively  over  the  period  of  the  binary  symbol,  m.  Bit 
decision  is  made  at  the  end  of  the  symbol  period.  As  in  [1],  decision- 
direction  is  to  be  used  from  symbol  to  symbol,  in  order  to  preclude  a 
processor  size  which  would  grow  exponentially  with  symbol  sequence  length. 

The  assumed  data  generating  model  is  that  of  Figure  3,  wherein 
z(k),  n(k),  4(k;m),  and  ^(k)  are  the  sampled  versions  of  z(t),  n(t), 

4(t;m),  andy(t),  respectively,  and  k  is  sample  number.  The  colored 
interference  process,  y.(k),  is  generated  from  zero-mean,  white,  Gaussian, 
unit-variance  noise  (a  two-vector),  W(k),  which  is  independent  of  the 
channel  noise,  n(k).  The  true  structure  of  the  ,j/(k)  generator  is  the  set 
{r  ,<t>,A}  which  may  also  be  unknown.  The  problem  of  joint  identification 
of  {r,4>,A}  has  been  treated  in  Reference  1. 

The  decision  on  m  is  to  be  made  according  to  the  maximum  A  Posteriori 
Probability  (MAP)  strategy.  That  is,  a  decision  statistic,  S  is  to  be 
formed  recursively  from  the  jj(|t  of  all  data  samples,  z(k),  taken  in  sequence 
during  the  symbol  period.  Let  Z(k)  denote  the  2-K  vector  of  K  samples  of 
the  I-Q  data  during  the^riod. 

T 

Z(k)  =  [z(K),  z(K-l),...,  z(l)]' 


z(K) 

Z(K-1) 


(7) 
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Figure  3.  Data  Generating  Model 


The  MAP  decision  statistic  is  the  probability 

S1 (K,m)  =  p(m|Z(K))  (8) 

A 

The  decision  rule  is  that  the  detected  symbol,  m,  is  that  one  for  which 
S  (K,m)  is  maximum. 

Assuming  that  the  A  Priori  probability  of  transmitted  symbols,  p(m), 
is  known,  maximization  of  S^(K,m)  is  obtained  by  just  maximizing  the 
Maximum  Likelihood  (ML)  statistic,  S(K,m),  where 

sV,m)  =  •  p(Z(K)|m) 

S(K,m)  =  p(Z(K)|m)  (9) 

Now,  the  signal,  a(k;m),  is  a  function  of  an  unknown  phase  process, 

(j>  (k),  as  per  Eq.  (4).  Thus,  define  a  K-vector,  $(K),  as 


6 


The  unknown  phase  process,  £(K),  is  imbedded  in  the  problem  by  using  the 
composite  detection  approach,  as 


p(Z(K)|m)  =  //•••/  p(Z(K),  i(K)|m)d<j>4(K)...d<|>s(l)  (11) 

The  ML  decision  statistic,  S(K,m)  is  to  be  generated  in  recursive  form. 
Thus,  the  argument  of  the  integral  in  (11)  is  manipulated  to  obtain  a 
recursive  form. 

We  have 

p(Z(K),  ${K)|m)  = 

=  p(z(K),  Z(K-l),  <J>4(K),  *(K-l)|m) 

-  p(z(K),  4>a(K)  |Z(K-1 ),  *(K-1),  m)  • 
p(Z(K-l),  £(K-l)|m) 

=  p(z(K),  <^(K)|Z(K-1),  i(K-l),  m)  • 
p(i(K-l)|Z(K-l),  m)  •  p(Z(K-l)|m)  (12) 

Then, 

p(Z(K)|m) 

=  //•••/  p(z(K) ,  <b  (K)|Z(K-1),  i(K-l),  m)  • 

<o 

p(i(K-l)|Z(K-l),  m)  •  p(Z(K-l)|m)d^(K)--*d^(l) 

(13) 

and 

S(K,m)  =  S(K-1 ,  m)Q(K)  (14) 

Q(K)  =  //•••/  p(z(K)|^(K),  £(K-1),  Z(K-l),  m)  • 

p(^(K)|£(K-l),  Z(K-l),  m)  •  p(£(K-l)|Z(K-l),  m)  • 

d^(K).--d^(l)  (15) 


where 


Now,  let  us  define  $(£)  to  be  the  conditional -mean  estimate  of  4>(£), 
given  the  data  zjk)  for  k  =  1,2,...,£,  and  given  the  symbol,  m.  Then, 
4>U)  maximizes  p(${£)|Z(£),  m).  Now,  it  is  assumed  that  the  gradients 
of  p(z(K)|^(K),  i(K-l),  Z(K-l),  m)  and  of  p(^(K)  |*(K-1 ) ,  Z(K-l),  m), 
with  respect  to  <px(K-l),...,  4>  (1),  evaluated  in  the  neighborhood  of 
$(K-1),  are  sufficiently  small  so  that  the  approximation  may  be  made 


Q(K)  «  /p(z(K)|4>4(K),  i(K-l),  Z(K-l),  m)  ‘ 

P<4>a(K)|Z(K-1),  i(K-l),  m)d^(K)  (16) 

This  approximation  says  that  the  functions  p(z(K)|(  ))  and  p(4>  (K)  (  )), 
viewed  as  functions  of  the  4>.{K-1),  •••*  <J>,(1)»  are  sufficiently  "flat" 
that  p($(K-l)|(  ))  appears  as  a  multi-dimensional  delta  function,  cen- 

A  A 

tered  at  the  co-ordinates,  (K-l ),. . . ,  ( 1 ) .  The  multiple  integral  then 

simply  evaluates  the  argument  at  those  coordinates,  analagous  to  "sifting" 
with  a  delta  function. 

Physically,  the  approximation  means  the  following.  If  a  sufficiently 
accurate  conditional -mean  estimate  may  be  obtained  for  the  phase  process, 

♦  .  0 ) t  ■  •  • , .  ( K— 1 ) ,  then  the  density,  p(4>(K-l)|Z(K-l),m),  will  have  a 
very  small  variance  about  the  mean  estimate.  Thus,  the  density  p(4>(K-l)| 

(  ))  will  be  so  highly  concentrated  that  the  densities,  p(z_(K)|(  ))  and 
p ( ( K )  |  (  ))  will  be  flat  by  comparison.  Thus,  the  accuracy  of  the  ap¬ 
proximation  depends  entirely  on  the  availability  of  a  very  good  phase 
estimate. 

A 

Similarly,  now  define  <p.(£)  to  be  the  one-stage  conditional -mean 
prediction  of  <j>.(K),  given  the  previous  data,  Z(£-l),  the  previous  condi- 

•6  A 

tional  mean  estimate,  $(£-1),  and  the  symbol,  m.  As  previously,  assume 

A 

that  4>  It)  is  sufficiently  accurate  so  that  p(z(£)|(  ))  is  flat,  by  compar- 

A  A 

ison,  in  the  neighborhood  of  <j>,{£).  This,  then,  yields  the  final  approxi- 
mation 

Q(K)  -  p(z(K)|J4(K),  i(K-l),  Z(K-l),  m)  (17) 

The  recursive  decision  statistic  is  then 


8 


08) 


K 

S(K,m)  =  tt  Q(k) 
k=l 


K  „ 

=  TT  p(z(k)|<t>  (k),  $(k-l),  Z(k-l),  m) 
k=l  “  4 

It  is  seen  from  (17)  and  (18)  that  the  recursive  detector  must  form  the 

A  A 

conditional  probability  function,  p(z(k)  |<f>  (k) ,  4>(k-l),  Z(k-l),  m),  at 
each  sample  time  (number)  k.  Moreover,  operating  in  parallel  with  the 
decision  circuitry,  and  furnishing  recursive  phase  estimates  to  it,  is 
a  conditional-mean  phase  estimator-predictor.  The  estimator  produces  the 
estimates 

^(k)  -  E{J4(k)|i(k-l),  Z(k-l),  m) 

i(k)  »  E{£(k)|Z(k),  m}  (19) 

The  problem  of  conditional -mean  estimation  of  the  phase  of  a  sinusoid 
in  Gaussian  noise  is  a  non-linear  estimation  problem  without  a  known 
general  solution.  However,  the  first-order  approximate  solution  is  known 
and  is  a  phase-locked  loop  [2].  The  closely  related  approximate  Maximum 
A  Posteriori  Probability  estimator  is  also  a  phase-locked  loop  [3].  Given 
the  symbol,  m,  and,  hence,  the  corresponding  signal  waveform,  *(t;m), 
the  bandpass  received  data,  z(t),  consists  of  a  sine  wave  of  unknown 
(random)  phase,  imbedded  in  additive  colored  plus  white  Gaussian  noise. 

Thus,  the  available  solution  t."  the  estimation  problem  indicated  by  (14) 
is  the  decision-directed  phase-locked  loop.  Note  that  the  PLL  is  only 
the  approximate  solution  to  (19)  for  the  case  where  the  phase-estimation 
error  is  quite  small.  Thus,  the  optimality  of  the  detection  algorithm  of 
(18)  will  depend  on  the  phase  estimation  accuracy  which  may  be  realized 
in  practice  using  the  PLL. 

3.  THE  I-Q  DATA  MODEL  WITH  PHASE  ESTIMATION 

In  order  to  proceed  with  the  detection  and  phase  estimation  algorithms, 
the  discrete-time  I-Q  data  generation  model  must  be  extended  beyond  that 
of  equation  (6)  and  Figure  3.  Under  the  assumption  that  the  I-Q  demodu¬ 
lating  reference  sinusoid  phases  are  estimated,  the  model  changes  somewhat. 
Let  the  physical  model  be  shown  in  Figure  4. 
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Figure  4.  Data  Model 


In  Figure  4,  the  transmitted  signal  with  unknown  phase  is 
4'(t,m)  =  Acos[u)  t  +  (fi(t,m)  +  <j>  (t)]  (20) 

C  <6 

where  <p{ t ,m)  is  the  angle  modulation  waveform,  containing  the  symbol,  m. 
The  unknown,  possibly  time-varying,  phase  term  is  <|>.(t).  The  additive, 
zero-mean,  Gaussian  colored  and  white  noises  are  respectively, 

y'(t)  =  y^(t)cos<oct  -  yq(t)sincjct 

n'(t)  =  nl(t)cosu)ct  -  n^(t)sinwct  (21) 

where  the  i  and  q  subscripts  denote  "in-phase"  and  "quadrature"  low-pass 
components,  respectively. 

The  product  detector  reference  sinusoids  are 

*ri(t)  =  2cos[wct  +  <{>4(t)] 

4rq(*)  =  -2sin[<i)ct  +  4>4(t)]  (22) 
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where  4>  (t)  is  the  phase  estimate  of  <j>  ft),  provided  by  the  phase-locked 
loop.  The  usual  problem  of  the  phase-locked  loop  responding  to  the  low 
frequency  portion  of  the  modulation  <f>(t,m)  may  be  encountered,  depending 
on  the  exact  form  of  the  modulation. 

Now  define. 


ijt)  -  <pjt)  &  e(t) 


(23) 


It  may  be  shown  that  the  low-pass  I-Q  data  vector  has  the  form 


Z^t) 

r 

_2q(t)_ 

L 

where 


cosE(t)  sine(t)-!  Acos4>(t,m)~j 

-sine(t)  cose(t)J  Asin<(>(t,m) 


yi(tri 

yq(t) 


n.(t)' 

nq(t> 


(24) 


y.,-  ( t) 

"cos 4>  (t)  sin^Jt)-" 

-6  4 

y\(t) 

yq(t) 

_-sind)  (t)  coso.U)- 

4  4 

”n.(tf 

rcos<f),(t)  sin<{>  (t)~ 

4  4 

"n:(t)_ 

nq(t)J 

_-sindi,(t)  cos<)),(t)_ 

4  -A 

(25) 


With  n'(t)  white,  Gaussian,  zero-mean  with  variance,  of,  then  n(t)  is 

2  r 

also  white,  zero-mean,  with  variance  ar-  This  is  because  the  multiplying 
matrix  is  a  rotation  matrix.  However,  njt)  is  not  Gaussian,  in  general. 
For  time  periods  which  are  short  compared  to  the  reciprocal  bandwidth 

A 

of  4>.(t),  n(t)  appears  approximately  Gaussian.  With^'(t)  colored, 

^  p 

Gaussian,  zero-mean,  with  variance  a  ,  ^(t)  is  zero-mean  with  variance 
2  * 

°y-  lit)  is  not  Gaussian  and  may  be  of  slightly  greater  bandwidth  than 

A 

j£'(t),  if  the  variation  of  4>4(t)  is  not  small. 

The  new  data  model  of  (24)  may  be  written  in  three  equivalent  forms, 
and  in  discrete  time,  as 
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z(  k)  =  H[^(k)lH[^(k)l»(k;m)  +  *'(k)  +  n’(k)}  (26a) 
z(k)  =  H[e(k)]4(k;m)  +  y(k)  +  n(k)  (26b) 
z.(k)  =  H[k;m]£(k)  +  ^(k)  +  n(k)  (26c) 


where  in  (26) 


H[e(k)]  = 


cose(k)  sinc(k)- 
L-sin£(k)  cose(k)_ 


H[k;m]  = 


cos<j>(k;m)  -sin<|)(k;m) 
Lsin4>(k;m)  cos<j>(k 


"’l, 

;m)  J 


A(k;m)  = 

£(k)  = 


;m)~ 

;m)_ 


Acos<j>(k;m) 
LAsin<j>(k 

~  Acose(k)“| 
_-Asine(k) J 


rcos*  (k) 
H  [  <0  ( k )  3  =  4 

Lsin<j>4(k) 


-sln^(k) 
cos<j>4(k) 
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(27) 


H[04(k)]  = 


/N 

r  cos<j>4(k) 

L-si 


1  ntj>4  ( k ) 


sin4>4(k)- 

/v 

cosi{>4(k)_ 


In  (26b),  the  matrix  H(k;m)  is  a  function  only  of  the  signal.  The  vector, 
£(k),  is  a  function  only  of  the  phase-tracking  error  process,  e(k) . 
Detection  of  m  in  the  presence  of  j>(k)  is  a  multiplicative  noise  detection 
problem.  The  presence  of  the  additive  colored  and  white  noise  processes, 
jr(k)  and  £(k),  respectively,  gives  a  compound  detection  problem,  having 
multiplicative  and  additive  colored  noise. 

The  compound  detection  problem  for  multiplicative  and  additive 
colored  Gaussian  noise  was  solved  in  [4].  There  it  was  found  that  the 
detector  was  one  which  tracked  both  the  multiplicative  and  additive  color¬ 
ed  noises  and  attempted  to  remove  them  from  the  data,  zU).  Although,  in 
the  present  case,  the  various  multiplicative  and  additive  noises  are  not 
strictly  Gaussian,  the  tracking  detector  may  still  be  used.  Note  that 
when  e(k)  is  small  then  £(k)  is  approximately 


p(k)  =  A[_e|kj]  ;  |e(k) j  «  1  (£8) 

In  this  case,  e(k),  the  phase  tracking  error,  is  Gaussian  and  £(k)  is 
approximately  Gaussian. 

The  final  data  generator  diagram,  corresponding  to  equations  (26) 
is  shown  in  Figure  5. 
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Figure  5.  Data  Generator  Model  for  Phase  Estimation 


4.  DETECTOR  STRUCTURE  AND  ALGORITHMS 

With  the  data  generator  given  as  in  Figure  5,  the  tracking  detector, 
with  phase  estimator,  takes  the  form  of  Figure  6.  In  the  detector,  there 
are  two  decision -directed  tracking  filters,  one  implemented  for  the 
signal  waveform  corresponding  to  m=0,  and  the  other  for  m=l .  Each  tracking 
filter  is  matched,  in  the  Wiener  sense,  to  both  £(k),  the  multiplicative 
noise,  and  y(k),  the  additive  noise.  Thus,  the  detectors  are  implemented 
for  the  data,  z_( k ) ,  in  the  form  of  equation  (26c).  The  tracking  error 
waveforms,  jj(k;m),  drive  the  decision  circuitry  which  produces  the  deci¬ 
sion  on  the  received  symbol  as  m. 

It  was  shown  above  that  generally  the  phase  estimator  is  decision- 
directed.  However,  a  non-decision-directed  phase  estimator  may  be  imple¬ 
mented  if  the  transmitted  signal  possesses  a  residual  unmodulated  carrier 
component.  This  is  shown  as  follows  for  a  phase-shift-keyed  signal. 
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Decision  Direction 


Figure  6.  Compound  Detector  and  Phase  Estimation 


Suppose  that  the  signal  phase  term  is 


Then 


<j>{k;m)  =  A<j>*c(k;m)  ;  c(k;m)  =  1;  m=0 

=  -1;  m=l 

0  <  A<j>  <  tt/2 


cosg>(k;m)  =  cos(A0) 
sin<j>(k;m)  =  c(k;m)  •  sin(A4>) 


(29) 


(30) 


It  follows  that 


H[k;m]^(k)  =  Acos(A4>) 


”  cose(k)~ 
_-sine(k)_ 


+ 


c(k;m)*Asin(A<}>) 


sine(k)“ 

_cose(k)_ 


(31) 


From  (31)  it  is  seen  that  there  is  present  in  the  received  data  an  addi¬ 
tive  term  proportional  to  -sine(k),  which  may  be  used  to  drive  the  phase 
estimator.  Likewise,  there  is  an  additive  term  proportional  to  cose(k) 
which  may  be  used  to  estimate  A  (coherent  automatic  gain  control).  The 
PSK  waveform,  c(k;m),  is  present  in  both  I-Q  channels,  due  to  the  multi¬ 
plicative  process  with  components  sine(k)  and  cose(k).  Provided  that  the 
bandwidth  of  c(k;m)  is  sufficiently  wide  and  the  closed-loop  tracking 
bandwidth  of  the  phase  estimator  is  sufficiently  small,  the  estimator  can 

track  phase  in  the  presence  of  c(k;m)  without  decision-direction. 

Each  decision-directed  tracking  filter  in  Figure  6  is  of  the  form 

of  Figure  7.  In  the  figure,  the  inner  loop,  composed  of  elements  G  , 

4>p,  Ap,  and  H[k;m],  track  the  multiplicative  process,  £(k).  The  elements 
(Gp,  4>p,  Ap}  are  the  elements  of  a  Wiener  filter  in  Kalman  canonical  form, 
matched  to  |>(k).  H[k,m]  contains  the  signal  waveform  elements,  as  in  (27). 
The  outer  loop  tracks  the  additive  colored  interference,  ^(k).  The  ele¬ 
ments,  {Gj,  Aj},  are  those  of  a  Wiener  filter  matched  to  y(k).  The 
filter  algorithms  are 


15 


Figure  7.  Tracking  Filter 
£(k)  =  ApXp(kjk-l ) 

Xp(k|k-1)  =  $p[Xp(k-l  | k-2)  +  G^k-l)] 
i(k)  =  Aj£j(k|k-1 ) 

£j(k|k-l )  =  ^j[Xj(k-l  | k-2)  +  G^k-l)] 

£(k)  =  z(k)  -  [H(k;m)£(k)  +  j(.(k)]  (32) 

A  A 

It  is  seen  from  Figure  7  and  (32)  that  the  ^(k)  filter  and  £(k)  filter 
are  uncoupled,  except  for  that  coupling  inherent  in  the  pseudo-innovations, 
£(k).  Filter  design  consists  of  selecting  the  two  sets  of  parameters 
(G.,  <&.,  A.)  and  {G  ,  $D,  A0}.  The  selection  is  based  on  either  real-time 
identification  of  ^(k)  and  £(k),  as  per  [1].  or  on  an  ad  hoc  worst  case 
design.  The  ad  hoc  design,  while  not  optimum,  would,  under  conditions 
discussed  in  [1],  produce  acceptable  results. 
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5.  THE  PHASE  ESTIMATOR 

From  equations  (26c)  and  (31)  we  may  write  an  expression  for  the 
(continuous-time)  data  vector,  as  seen  by  the  phase  estimator,  as 

-cose(t)- 

z(t)  =  A'  +  n(t)  (33) 

_-sine(t)_ 

In  (33),  ji( t )  is  the  total  noise  process  due  to  ^(t),  njt),  and  c(t:m). 

For  the  bandwidth  of  y(t)  and  the  band-rate  of  c(t;m)  sufficiently  great 
with  respect  to  the  closed  loop  bandwidth  of  the  phase  estimator,  the 
noise  process,  n(t),  will  appear  white  to  the  phase  estimator. 

A 

It  is  seen  that  the  problem  of  deriving  the  phase  reference,  <p  (t), 

•6 

which  is  an  accurate  estimate  of  the  residual  carrier  phase,  <j>  (t),  is 
that  of  minimising  e(t)  in  the  presence  of  the  unknown  amplitude.  A',  and 
noise,  n(t).  This  is,  essentially,  a  phase-locked  loop  problem.  Under 
the  assumption  that  n(t)  is  white  and  Gaussian,  the  solution  is  the  clas¬ 
sical  phase-locked  loop. 

Note  that  the  usual  problem  of  unknown  signal  amplitude.  A’,  is 
present.  There  are  two  classical  solutions.  One  is  to  use  the  Q-channel 
only,  for  phase  estimation,  with  an  ideal  pre-limiter  to  remove  dependence 
on  A'.  The  other  solution  is  to  also  use  the  I-channel  to  estimate  A' 
and  to  then  control  the  gain  of  the  Q-channel.  An  extension  of  the  second 
method  is  shown  in  Figure  8. 

In  Figure  8,  the  Q-channel  waveform,  zq(t)  is  processed  by  a  "Loop 
Filter"  with  low  frequency  gain,  H(0),  to  produce  an  estimate  of  the  term, 
(-A'sine(t)),  weighted  by  H(0).  The  I-channel  waveform,  (t) ,  is  process 
ed  by  a  low-pass  filter  with  unit  low  frequency  gain  to  procuce  an  esti¬ 
mate  of  the  term,  A'cose(t).  The  two  filter  output  terms  are  then  divided 
point-wise  in  a  digital  divider  to  provide  an  estimate  of  (-tane(t)), 
weighted  by  H(0).  The  latter  estimate  then  drives  the  Voltage-Controlled- 

A 

Oscillator  (VCO)  to  produce  the  reference  phase,  <f\(t).  It  can  be  seen 

•6 

from  the  defining  equation  (23)  for  e(t)  that  the  mechanization  of 

A 

Figure  8  causes  4>.(t)  to  track  <j>  (t). 


Zq(f ) 


Zj(t) 


Figure  8.  Tangent  Phase-Locked  Loop 


The  usual  phase-locked  loop  generates  a  tracking  error  voltage  pro¬ 
portional  to  (-sine(t)).  The  present  implementation  provides  a  tracking 
error  proportional  to  (-tane(t)),  which  will  yield  higher  loop  gain  for  a 
large  tracking  error,  e(t).  However,  the  main  reason  for  using  the  "Tan¬ 
gent-Loop"  mechanization  is  to  obtain  the  automatic  gain  control  feature 
in  the  cancellation  of  the  unknown  amplitude.  A'. 

The  design  of  the  loop  parameters,  notably  the  loop  filter,  is  per¬ 
formed  by  assuming  linear  operation  of  the  loop.  That  is,  when  e(t) 
is  small,  say  less  than  12°  in  magnitude,  then  the  approximation  holds 

tane(t)  =  sine(t)  =  e(t)  (34) 

Then,  the  overall  system  operates  as  a  linear  servo-mechanism  for  phase, 
or  as  a  linear  phase-locked  loop. 

In  the  usual  implementation,  the  Loop  Filter  in  the  quadrature  channel 
is  implemented  with  one  finite  zero  of  transmission  and  one  finite,  non¬ 
zero,  pole.  The  pole  frequency,  zero  frequency,  and  low-frequency  gain, 

H ( 0 ) ,  are  set  to  realize  the  desired  closed-loop  noise  bandwidth,  static 
phase  error  for  VCO  frequency  offset,  and  second  order  dynamic  response. 

The  low  pass  filter  in  the  I-channel  is  set  for  the  same  zero  and  pole 

• 

frequencies  as  for  the  Q-channel  Loop  Filter,  but  with  unit  low-frequency 
gain. 

Note  that  for  the  PLL  to  operate  properly,  the  signal  to  noise  ratio 
must  be  large  in  the  closed-loop  equivalent  noise  bandwidth  of  the  loop, 
itself  The  PLL  bandwidth  is  to  be  maintained  small  enough  to  just 
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accomodate  the  dynamics  of  the  received  signal  phase,  $4(t),  due  to  Doppler 
effects  on  the  transmission  link.  For  the  case  where  the  incident  noise  is 
dominated  by  colored  interference,  such  as  jamming,  the  loop  perfr  •, nance 
will  be  affected  by  that  portion  of  the  colored  interference  falling  with¬ 
in  the  (narrow)  loop  bandwidth. 

6.  THE  LOOP  FILTER  MECHANIZATION 

The  continuous-time  version  of  the  Loop  Filter  is  characterized  by 
the  transfer  function 

H(s)  =  Kt-g]  (35) 

where  K,  z,  and  p  are  real,  with  z  and  p  being  negative.  Let  y(t)  and 
z(t)  denote  the  filter  input  and  output,  respectively.  A  state  variable 
representation  is  set  up,  using  the  single  filter  state,  x(t),  as 

x(t)  =  px ( t )  +  y(t) 

z(t)  =  K(p-z)x(t)  +  Ky(t)  (36) 

The  filter  is  converted  to  discrete  time  by  driving  it  with  an  ideal 
sampler  and  zero-order  hold  circuit  and  observing  the  output  only  at 

sampling  instants,  t  =  tk  for  k  =  1,2,3 .  The  differential  equation  of 

(36)  is  then  solved  between  the  kth  and  (k+l)st  sampling  times  as 

x((k+l)T)  =  exp[p( (k+1 )T-kT)]  •  x(kT) 

+  kT/(k+1)T  exp[p(k+l )T  -  T]W(T)dT  (37) 

where 

W(t)  =  p(kT);  kT  <  t  <  (k+l)T  (38) 

and  T  is  the  sampling  interval.  The  differential  equation  solution  then 
yeilds  the  governing  difference  equation  (discrete-time)  for  the  filter  as 
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f 


x(k+l)  =  <j>x(k)  +  yp(k) 
z(k)  =  K{p-z)x(k)  +  Kp(k) 

<t>  =  exp(pT)  :  y  =  l/p(<H)  (39) 


The  Loop  Filter  constants,  K,  z,  p,  are  set  according  to  specifica¬ 
tions  on  the  linearized  closed-loop  transfer  function  for  phase.  The  VCO 

A 

output  phase,  <fs,(t)  is  given  by 
& 


or 


<f^(t)  =  /{-[(j>4(t)  -  <j>4(t)]*h(t)}dt 

O .U)  -  *.U)]-H(s) 

A  t  ~  4 _ 4 _ 


(40) 


(41) 


where  H(s)  is  the  Loop  Filter  transfer  function  given  in  (35). 
The  closed-loop  transfer  function  for  the  PLL  is  then 


G(s) 


=  H(s) 

^(s)  s  + H(s) 


Substituting  for  H(s)  yields 


(42) 


G(s)  =  ,  K<w> -  =  ■,  Kts~z> - ,  (43) 

s  +  (K-p)s  -  Kz  s  +  26wns  +  </ 

where  6  and  are  the  classical  damping  ratio  and  resonant  frequency  for 
a  second-order  servo  system. 

The  Loop  Filter  low  frequency  gain,  H(0),  is  given  by 

H(°)  =  5^  K(|f|)  =  K  f  (44) 

For  most  PLL  designs  the  following  assumptions  hold 


-z  «  H (0 ) 

-p  «  K  (45) 

Thus,  by  equating  like  terms  in  the  denominator  of  (43) 
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(46) 


K  *  26w 


-Kz  =  u). 


Now,  it  may  be  shown  that  the  one-sided  closed- loop  noise  bandwidth, 
in  Hz,  for  G(s)  is  [5] 


Thus, 


p  _  KrK-z-i  K-z  _  wn  r,  .  . ,.2-. 

Bn  ~  4[K^J  -  T  ~  86  [1  +  46  ] 

r  2 

K  =  •  B 

_1  +  4  62J 

2  =  -  f - — ol  •  B 

_1  +  4  62J  n 


For  loop  dynamic  stability,  the  damping  ratio  is  set  as 


K  =  8/3  B 


z  =  -  4/3  Bn  =  -K/2  (50) 

The  Loop  Filter  pole  frequency,  p,  is  generally  set  as  small  as 
possible  in  magnitude.  This  is  because  p  affects  the  "static  phase  error" 
when  tracking  with  a  fixed  Doppler  offset  in  the  received  frequency.  In 
order  to  hold  the  loop  in  lock  when  the  input  phase  <{>  (t)  has  a  constant 
first  derivative  requires  a  constant  driving  voltage  into  the  VCO  and  hence 
a  constant  phase  error,  e(t).  Thus, 

~  i  (t)  ^  Aw  =  -  H(0)  tane$p  (51) 

where  esp  is  the  static  phase  error  for  a  Doppler  offset,  Au>  =  2nAf.  The 
d.c.  gain  of  the  loop  filter  is 
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(52) 


(53) 


(54) 


Equation  (54)  gives  the  relation  between  the  various  quantities  and  f  . 

Thus,  the  design  equations  for  the  quadrature  channel  loop  filter 
are 

K  =8/38 
q  n 

z  =  -4/3  ;  quadrature  filter  (55) 


where  esp  is  static  phase  error  in  radians  for  a  Doppler  offset  of  Af  Hertz 

and  a  closed  loop  noise  bandwidth  of  B  Hz. 

n 

For  the  inphase  filter,  the  same  pole,  p,  and  zero,  z,  are  used,  but  the 
d.c.  gain  is  reduced  to  unity  to  give  a  filter  gain  constant 

K.  =  p/z  :  inphase  filter  (56) 

The  block  diagram  of  the  phase  estimator  is  given  in  Figure  9.  In 
the  figure,  the  discrete  time  version  of  the  VCO  (phase  integrator)  is 
represented  by 

JA(k+l)  =  J4(k)  +  T/2[v(k+l )  +  v(k)]  (57) 

where  v(k)  is  the  VCO  input. 


22 


SECTION  III 


ON  THE  EXISTANCE  OF  NON-COHERENT  TRACKING  DETECTORS 

It  is  desired  now  to  determine  if  a  non-coherent  version  of  the 
tracking  detector  exists.  In  [1]  the  non-coherent  version  of  the  standard 
FSK  detector  for  white  noise  was  derived.  The  approach  for  the  tracking 
detector  will  be  similar.  An  unknown  constant  phase  term  will  be  intro¬ 
duced  into  the  formulation  of  the  detection  problem.  Then,  the  detec¬ 
tion  statistic  will  be  averaged  with  respect  to  the  unknown  phase.  Up 
to  this  point,  the  procedure  is  the  same  as  was  followed  in  1 1. 2.  That  is, 
the  problem  is  that  of  composite  detection  for  unknown  phase.  In  II. 2 
there  existed  a  solution  of  the  composite  detection  problem  which  produced 
a  phase  estimator  as  part  of  the  detector.  In  the  present  formulation, 
the  phase  estimator  solution  is  purposely  rejected  and  no  attempt  is 
made  to  take  advantage  of  possible  phase  information.  Rather  the  unknown 
phase  is  defined  to  be  uniformly  distributed  over  the  interval,  [0,  2tt], 
and  to  be  a  constant  random  variable  over  the  time  interval  of  the  signal 
symbol.  Then  it  is  to  be  determined  whether  averaging  the  decision  statis¬ 
tic  over  phase  produces  a  sufficient  statistic  for  detection. 

The  unknown  phase  enters  the  problem  as  per  Figure  2,  where  now 
cf>o(t)  is  defined  to  be  constant  over  the  symbol  interval,  which  is  also 
the  processing  time.  Also  <t>Q(t)  is  uniformly  distributed  as 


4>0(t)  =  4> 

p((t>)  =  1/2tt 
=  0 


:  0  <  t  <  T 

:  0  <_<)><_  2n 

otherwi se 


(58) 


The  discrete  time  data  model,  z(k),  is  essentially  that  of  (26a) 

A 

where  <j>  (k)  =  <p  and  <j>  (k)  =  0.  Thus, 

c  0  <6 

z(k)  =  H(4>)[<i(k)  +  y(k)  +  n(k)]  (59) 


where  ajk)  is  the  transmitted  signal,  /(k)  is  the  colored  interference, 
and  n(k)  is  the  white  noise. 

The  detection  statistic,  S(K),  is  formed  recursively  from  the  z(k), 
and  is  the  Maximum  A  Posteriori  Probability  function,  p(m|z(K)),  where 
Z(K)  is  the  2K  partitioned  vector. 
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(60) 


Z(K)  =  [z(K),  z(K-l),...,  z(l)]T 

The  quantity,  m,  is  the  signal  digit,  which  for  the  binary  case  is  either 
0  or  1.  Under  the  assumption  that  the  transmitted  digits,  m,  are  equally 
distributed  (p(m)  =  1/2;  m=0,l),  the  MAP  statistic  is  equivalent  to  the 
Maximum-Likelihood  (ML)  statistic,  p(ZjK)|m).  Thus,  S(K)  is  obtained  by 
averaging  the  joint  density  on  Z(K)  and  <p,  given  m. 

2tt 

S(K)  =  p(Z(K)|m)  =  /  p(Z(K),  4>|m)d<|> 

0 

2n  , 

=  /  ^  P(Z(K)  |m,  4>)d<j) 

The  conditional  density,  p(Z(K)jm,  <)>)  is 
K 

p(Z(K)|m,4>)  =  -rr  p(z(k) JZ(k-l ) ,  m,  <j>) 
k=l 


(61) 


(62) 


Now,  p(z(k)  |Z(k-l ),  m,  <j>)  is  Gaussian,  under  the  definition  that  y(k)  and 
n(k)  are  Gaussian,  and  is  given  by 

p(z(k)  |Z(k-l),  m,  <ji)  = 


=  exp[- 

2-rra  2a 

v  v 


(z(k)  -  z(k|k-l,  m,  «p) )T(z^(k)  -  z(k|k-l,  m,  <(>))] 

(63) 


2  ^ 

In  (63),  is  the  steady-state  Innovations  variance  and  zjkjk-1,  m,  )  is 

the  recursive  estimate  of  the  kth  data  sample,  given  all  the  data  up 
through  the  (k-l)st  sample.  This  one-sample  predictive  estimate  is  ob¬ 
tained  from  the  Kalman-form  filter  of  Figure  10.  In  the  figure,  the 
quantities,  {$,  A,  G},  are  the  appropriate  Kalman  (Wiener)  filter  para¬ 
meters  for  tracking  y(k),  the  colored  interference,  in  the  presence  of 
n(k),  the  white  noise. 
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Figure  10.  Kalman  Filter 

The  filtering  algorithms  are 

z(k|k-l)  =  H(<t>)[a(k)  +  i(kjk-l)  =  H(<j>)a(k)  +  HU)a$xU-D 

x(k)  =  'fU)x(k-l)  +  n(k,<j>);  Tf(4>)  =  [I-GH(<j>)A]<i> 

y(M)  =  G[z(k)  -  H(tj>)a(k)]  (64) 

The  solution  to  (64)  at  the  kth  sample  is  given  by 

k-1 

z(k|k~l,  m,  <t>)  =  H(<J>)[a(k;m)  +  A4>[fk'1  (<j>)x(0)  +  l  (4>)y(k-i ) ,<)>)]] 

i=l 

(65) 

It  is  seen  at  this  point  that  any  hope  of  averaging  p(Z_(K) jm.tp)  over  <j>  is 

A 

futile  due  to  the  internal  dependency  of  z_(k j k-T ,  m,  <p)  on  .  That  is, 
it  is  the  feedback  dependency  of  the  estimate  y.(k)k-l,  m,  <t>)  upon  <}>  which 
defeats  the  prospect  of  averaging  over  4>. 
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s(k;m) 


Figure  11.  Kalman  Filter 

There  is  a  second  possibility  for  a  noncoherent  implementation.  The 
term,  H(<j>)^(k),  in  the  data  model  of  (59)  is  not  strictly  Gaussian,  but 
does  have  the  same  first  and  second  moments  as  yU),  since  H(<)>)  is  unitary. 
Also,  since  cf>  is  constant  over  a  symbol  period,  H(<j>)y(k)  has  the  same 
short-term  spectral  properties  as  y(k).  Thus,  the  data  form  may  be  re¬ 
defined  as 

zU)  =  H(<}>)4(k;m)  +  y^(k)  +  ji(k)  ;  1  £  k  £  K  (66) 

where  y(k)  and  ji(k)  have  replaced  H(4>)y(k)  and  H($)ji(k),  respectively.  In 
(66),  y(k)  and  n(k)  are  taken  as  Gaussian. 

The  resulting  Kalman  estimator  for  zjk|k-l,  m,  <f>),  corresponding  to 
the  data  model  of  (66)  is  as  in  Figure  11. 

The  filtering  algorithms  now  are 

A 

z(k|k-l,  m,  <p)  =  H(<j>)a(k;m)  +  y(k|k-l,  m,  <p) 

=  H(4>)4(k;m)  +  A$x(k-1) 

x_(k)  =  l'x(k-l)  +  jj( k,4>)  ;  V  =  (I-GA)4> 

u(k;$)  =  G[z(k)  -  H(<j>)a(k;m)]  (67) 
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The  solution  to  (67)  is 

z(k|k-l,  m,  <}>)  =  H(<j>)4(k;m)  +  A^O^xtO)  +  £  4,i _1\i(k-i  ,<{>)] 

i=l 

(68) 

Now,  (68)  is  somewhat  of  an  improvement  over  (65)  in  that  V  is  no  longer 
a  function  of  4>.  Unfortunately,  jj(  )  is  still  dependent  on  <j>  and  this 
causes  the  dependency  of  zjkjk-1,  m,  <f>)  on  $  to  be  internal  because  of 
the  feedback  structure  of  the  filter.  Thus,  averaging  p(_Z(K)  | m,g> )  over 
<J>  is  still  not  feasible. 

The  argument  of  the  exponent  of  p(z(k) |Z(K-1 ) ,  m,  <}>)  in  (63)  is 

Arg  =  (z(k)  -  y.(k|k-l,  m,  <p) )T( • )  +  4T(k;m)4(k;m) 

-  24T(k;m)HT(<|,)[z(k)  -  i(k|k-l,  m,  0)]  (69) 

This  argument  is  of  the  same  form  as  is  encountered  in  the  standard  non- 

A 

coherent  FSK  detector  problem  [l],  except  that  (z_(k)  -  ^(k|k-l,  m,  <b ) ) 
has  replaced  ,z(k).  Were  it  not  for  the  fact  that  x(k|k-l,  m,  $)  is  an 
explicit  function  of  <p,  as  in  (67),  then  the  averaging  over  4>  would  be 
exactly  the  same  as  in  the  FSK  problem.  Unfortunately,  there  seems  to 
be  no  further  recourse  to  the  problem  at  this  point. 
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SECTION  IV 


SIMULATION  RESULTS 

A  Monte-Carlo  simulation  program  was  written  to  obtain  error-rate 
results  for  coherent  detection  with  phase  estimation.  The  detecto  algo¬ 
rithm  which  was  implemented  was  that  detailed  in  Section  II. 4.  The  program 
realized  the  compound  detector  and  phase  estimator  of  Figure  6  where  the 
tracking  filters  were  of  the  form  given  in  Figure  7.  The  phase  estimator 
was  the  Tangent  Phase-locked  loop  shown  in  Figure  9. 

In  order  to  reduce  simulation  run  times,  the  Monte  Carlo  program, 
documented  in  [4],  was  not  modified  for  present  use.  Rather,  an  entirely 
new  program  was  written.  In  the  new  program,  the  data  generator,  shown 
in  Figure  5,  was  reduced  from  three  states,  as  in  [4],  to  one  state.  This 
resulted  in  the  Kalman  filters  also  having  one  state  in  each  branch 
shown  in  Figure  7.  Since  computation  load  increases  exponentially  with 
state  size,  a  considerable  savings  was  made.  All  that  was  lost  was  some 
flexibility  in  modeling  the  additive  colored  noise  process.  For  the  pur¬ 
poses  of  the  present  work,  the  one-state  model  was  sufficient. 

It  was  desired  to  test  the  compound  detector  and  phase  estimator  in 
a  realistic  but  stressful  environment.  Thus,  a  phase-locked  loop  noise 
bandwidth  of  2.5  ilz  was  chosen  as  being  as  small  as  could  likely  be  real¬ 
ized  in  a  reasonable  implementation.  It  was  desired  to  run  the  phase- 
locked  loop  at  0.3  radians  r.m.s.,  phase  error,  or  less.  Thus,  it  was 
necessary  to  relate  the  various  simulation  parameters,  such  as  E/NQ, 
colored  noise  bandwidth,  etc.,  to  the  phase-locked  loop  signal  to  noise 
ratio. 

Letting  J  denote  the  power  of  the  colored  process,  y(k),  (in  bandpass 
form)  and  Bj  the  one-sided  equivalent  noise  bandwidth  of  the  low-pass  I-Q 
process,  an  eqivalent  white  bandpass  spectral  density,  Nj,  for  the  colored 
process  is  defined  by 

J  =  Nj  •  2  Bj  (70) 

Then,  the  total  equivalnet  white  noise  spectral  density  is 


where  NQ  is  the  density  of  the  incident  additive  white  receiver  noise. 

The  symbol  energy,  E,  in  the  received  signal  is  related  to  total  sig¬ 
nal  power,  S,  and  symbol  period,  T,  by 

E  =  S*T*LM(A<f>)  (72) 

where  L^(A<|>)  is  the  "modulation  loss"  factor  given  by 

Lm(A*)  =  sin2(Ac|>)  (73) 

where  A<p  is  phase  deviation  in  radians  for  the  phase-shift  keyed  signal. 
Thus, 

C  _  E  /T. \ 


From  (70)  and  (74)  results 


<fo>"o s  E  ■  <!>VmM2  r>-No 


where  R  =  1/T  is  symbol  rate.  Thus, 
(E/NJ  n 

NJ  =  - V  •  (2T*)-No 

LM(Ad>)*(f)  2BJ  ° 


nt  =  [1  + 


(E/No)  R 

- V  •  (?in3N0 

Lm(A4>)-(§)  2BJ  0 


It  is  desired  to  compute  the  ratio  of  residual  carrier  power  to  total 
white  noise  power  in  the  Loop-noise  bandwidth  (one-sided),  B^.  The  resi¬ 
dual  carrier  power,  is 

L-(A4>)  £  nr 

-  Ec‘M>s  *  f ■  (79) 


flStwwS'*’ 
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where  L^(A4>)  is  “carrier  loss"  given  by 


L c  ( )  =  cos  (A<j>) 


The  desired  signal  to  noise  ratio  is 


( R/Bfj )  “  (  e/n  ) 


N  Bn  ~  NXBW  '  , 

'  N  T  N  tan2(A4>)  [1  + 


N;  '‘-'"o' 

^  (Xn 

- V-  •  ^2B/J 

L„(4»)-(f)  J 


Note  that  when  the  equivalent  white  spectral  density  of  the  colored  inter- 
ferrrrg  process  4s  much  larger  than  the  receiver  white  noise  spectral 
density,  then  (81)  reduces  to 


rjBN  ■ 2  m 

The  loop  phase  error  variance,  under  the  assumption  that  the  loop 
is  operating  linearly  for  phase  (large  loop  signal  to  noise  ratio),  is 
given  by 


2  __ 1 

°<P  "  Sr| 


radians 


and,  from  (83)  and  (81) 


2  2 
a  =  tan  (a$) 


1  2tVBJ*  ~l 

Lw 


Figure  12  shows  simulation  results  for  the  case  of  narrow-band  inter¬ 
ference  for  binary  phase-shift-keying  (PSK).  The  equivalent  square  band¬ 
width  of  the  colored  interference  process  is  275  HZ.  The  signal  symbol 
rate  is  2500  baud.  Thus  the  "bandwidth  to  bit-rate  ratio"  is  BW/BR  = 
0.109.  This  is  the  same  case  for  which  extensive  previous  results  were 
reported . 


Standard  Detector  - 
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with  Phase-Locked 
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Coherent  Detection  | 
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Perfect  Phase 
—  Reference 


SIMULATED  ERROR  RATE  -  COLORED  INTERFERENCE 

2500  Symbols/Second  10  Samples  /Symbol 
Binary  PSK,  Mod.  Index =  0.785 
Interference  Bandwidth  *  275  Hz. :  BW/BR*  0.109 
Multiplicative  Noise  Tracking  Bandwidth  «2.5 Hz. 

Phase -Locked  Loop  Bandwidth  =  2.5  Hz. 

Phase  Jitter  s  5  degrees,  r.m.s. 

Perfect  Identification 


SYMBOL  ENERGY  TO  NOISE  SPECTRAL  DENSITY  RATIO ,  E/N,,dB. 


Figure  12.  Simulation  Results 


For  the  case  of  Figure  12,  the  ratio  of  signal  power  to  additive 
colored  noise  is  unity,  or  zero  dB.  A  phase-locked  loop  is  implemented, 
as  described  in  Sections  II. 5.  and  I I. 6.  The  loop  noise  bandwidth  is  2.5 
Hz,  being  the  smallest  assumed  to  be  practical  for  this  case.  The  detec¬ 
tor  employs  both  additive  noise  tracking  and  multiplicative  noise  tracking 
with  the  latter  matched  to  a  multiplicative  noise  bandwidth  of  2.5  Hz. 
Perfect  identification  is  assumed  for  the  colored  additive  noise. 

From  (84),  the  predicted  value  of  loop  phase  jitter  is  determined  to 
be  5.4°  r.m.s.  The  actual  r.m.s.  values  recorded  in  the  simulation  were 
between  1.7°  and  9.6°  for  runs  up  to  1500  symbols  in  length.  The  loop  was 
observed  to  always  be  in  lock,  slipping  no  cycles  during  any  run. 

The  results  plotted  in  Figure  12  include  the  reference  graphs  of 
coherent  PSK  detection  for  white  noise  only,  and  IDEI  detection  with  per¬ 
fect  phase  reference.  Also,  is  shown  the  behavior  of  the  standard 
discrete-time  matched  filter  detector.  The  matched  filter  is  seen  to 
saturate  at  an  error  rate  of  0.14,  as  usual  [1],  The  IDEI  detector  is 
seen  to  yield  a  convex  error  rate  curve  for  -10  dB  <  E/Nq  £  20  dB.  How¬ 
ever,  for  20  dB  <  E/Nq,  the  slope  of  the  error  rate  curve  becomes  much 
less  steep.  Although  the  error-rate  continues  to  decrease  for  increasing 
E/Nq,  the  rate  of  decrease  is  not  as  good  as  was  obtained  for  ‘'pure" 
multiplicative  noise  in  [4], 

The  implications  (or  "cause")  of  the  change  in  slope  of  the  error  rate 
curve  for  20  dB  <  E/Nq  are,  at  present,  unknown.  Clearly,  there  is  a 
transition  at  E/Nq  -  20  dB  for  the  case  shown.  It  has  been  observed  in 
the  past  that  such  transitions  may  be  due  to  the  breakdown  of  basic 
modeling  assumptions  on  which  the  "optimum"  detector  is  founded.  One  such 
questionable  assumption  which  is  suspect  here  is  that  the  multiplicative 
noise  process,  due  to  carrier-tracking  phase  error,  is  Gaussian.  Also,  it 
may  be  that  the  phase-tracking  detection  algorithm  is  subject  to  an  ir¬ 
reducible  error-rate,  as  detailed  in  [8].  It  is  noted  that  the  IDEI  detec¬ 
tor  for  multiplicative  noise  has  not  previously  shown  such  an  irreducible 
error. 

In  conclusion,  this  simulation  for  the  SJR  =  0  dB  case  shows  that  much 
of  the  performance  measured  previously  for  perfect  phase  is  lost,  when  a 
standard  phase-locked  loop  is  used  in  parallel  with  the  IDEI  detector.  It 
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is  recalled  that  this  implementation  is  not  the  true  optimum,  for  two 
reasons.  One  is  the  Gaussian  multiplicative  noise  approximation.  The 
second  is  that  the  phase-tracking  loop  is  external  to  the  detector  and, 
thus,  does  not  take  advantage  of  the  colored  noise  tracking  capability  of 
the  detector  itself.  It  may  well  be  that  a  more  optimum  implementation 
will  result  by  imbedding  the  phase-estimation  algorithm  within  the  detector 
itself. 


SECTION  V 


COMPLETE  RECEIVER  ALGORITHMS 

1.  A  PROPOSED  BIT  SYNCHRONIZATION  ALGORITHM 

So  far  in  the  investigation  of  IDEI  detection,  it  has  been  assumed 
that  bit  timing  information  is  available.  This  is  important  for  the 
detector  in  terms  of  setting  the  start  and  stop  times  of  the  computation 
which  produces  the  decision  statistic,  S(K).  However,  now  the  synchro¬ 
nization  problem  is  finally  examined. 

Many  practical  bit  synchronizers  are  based  on  the  "Delay-lock  Loop," 
[6,  7].  This  technique  applies  to  any  coherent  signalling  scheme,  but 
is  generally  used  for  phase-shift-keying  (PSK).  Generally,  the  implemen¬ 
tation  uses  two  signal  cross-correlators  driven  with  time-staggered  signal 
reference  waveforms.  The  correlator  outputs  are  time- staggered  versions 
of  the  noisy  signal  autocorrelation  function.  By  subtracting  the  stag¬ 
gered  autocorrelation  functions,  a  tracking  error  function  is  produced 
which  drives  the  reference  genereator  into  bit  synchronism  with  the 
receiver  signal. 

The  key  to  the  functioning  of  the  del  ay- lock  bit  synchronizer  is 
the  production  of  a  signal  (from  the  correlator  output)  which  is  a  pos¬ 
itive,  even  function  whose  maximum  occurs  when  the  reference  generator  is 
in  synchronization  with  the  received  bit.  Those  positive  even  functions 
(autocorrelation  functions)  also  happen  to  be  the  sufficient  statistics 
for  detection  for  the  standard  detectors  which  use  delay-lack  bit  synchro¬ 
nization. 

In  the  IDEI  detector,  the  sufficient  statistic  for  detection  is  the 
pseudo-innovations  process,  or  noise  tracking  error.  It  was  seen  in  [1] 
that  there  was  associated  with  the  statistic  a  function  which  was  positive, 
with  minimum  value  occuring  for  perfect  identification  of  the  required 
noise  statistics.  With  "positive"  or  "negative"  identification  errors 
(in  the  sense  of  Figures  36  and  43  of  [1]),  the  function  value  increased. 
The  function  was  the  variance  of  the  noise  tracking  error. 

Now,  it  is  conjectured  that  the  IDEI  tracking  error  variance,  which 
is  necessarily  positive,  will  be  minimum  for  the  reference  signal,  a ( k ; n ) , 
exactly  synchronized  with  the  received  bit.  It  is  also  conjectured  that 
the  variance  will  increase  as  the  reference,  4(k;n),  becomes  unsynchro- 
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nized,  reqardless  of  whether  a(k;  n)  leads  or  lags  the  received  bit.  If 
this  conjecture  proves  true,  then  it  is  a  simple  matter  to  use  the  recip¬ 
rocal  of  the  tracking  error  variance  in  the  same  fashion  that  the  Delay- 
Lock  Loop  uses  the  autocorrelation  function,  to  form  a  synchronization 
tracking  error  function. 

2.  THE  COMPLETE  ALGORITHM 

The  complete  IDEI  algorithm  (excluding  identification)  can  be  postu¬ 
lated  as  follows,  for  binary  signalling.  See  Figure  13.  Two  IDEI  detec¬ 
tors,  with  imbedded  phase  estimators  are  implemented,  one  with  early 
waveform  reference  signals  and  one  with  late.  Each  detector  contains  two 
tracking  filters  of  the  form  of  Figure  7.  In  each  detector  are  produced 
the  detection  statistics,  SQ  and  S. ,  which  are  the  tracking  error 
variances,  conditioned  on  the  two  different  received  symbols,  m=0  and 
m=l,  respectively.  In  each  detector, symbol  decision  is  made  as  usual. 
Based  on  the  symbol  decision,  the  assumed  correct  tracking  error  variances, 
and  §£,  are  produced  by  the  early  and  late  detectors,  respectively. 

The  reciprocal  of  each  variance  is  taken  and  the  results  subtracted  to 
form  a  "Synch  Control"  driving  signal,  which  is  filtered  with  suitable 
gain  and  time  constant.  A  modulo-2  adder  is  implememted  to  determine  if 
the  decisions  in  the  early  and  late  detectors  do  not  result  in  the  same 
detected  symbol.  If  not,  the  synch,  control  signal  is  inhibited,  and 
synch  is  maintained  as  previously.  Decision-directed  reinitialization 
of  the  filters  is  carried  out  in  the  usual  manner,  independently  in  the 
early  and  late  detectors. 
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Synch.  Control  Signal 


Figure  13.  Complete  Algorithm  Oiogram. 


SECTION  VI 


CONCLUSIONS 

The  research  documented  in  this  report  has  yielded  several  interest¬ 
ing  results.  These  are  summarized  below  in  the  order  of  the  governing 
tasks  in  the  Contract  Statement  of  Work. 

Task  4. 

The  IDEI  (interference-tracking)  detection  algorithms  were  extended 
to  include  provision  of  the  required  carrier  phase  reference  through 
phase  tracking.  A  separate  phase-locked  Loop  was  implemented,  process¬ 
ing  the  received  data  in  parallel  with  the  detection  algorithm,  itself. 
The  detection  algorithm  was  augmented  to  track  the  multiplicative  noise 
resulting  from  the  phase  reference  variations,  as  well  as  tracking  the 
colored  additive  noise. 

Task  5. 

It  was  shown  analytically  that  a  non-coherent  version  of  the  IDEI 
detection  algorithm  does  not  exist.  This  result  is  due  to  the  feed¬ 
back  structure  inherent  in  the  IDEI  tracking  filter.  The  internal 
dependency  of  the  detection  statistic  on  the  unknown  phase  makes  it 
impractical  to  carry  out  the  phase  averaging  necessary  to  obtain  a  non¬ 
coherent  type  algorithm. 

Task  6. 

Based  on  the  result  of  Task  5,  a  non-coherent  IDEI  detector  for 
Differential  Phase  Shift  Keying  is  also  impractical  of  derivation. 

Task  7. 

A  bit  synchronization  technique  was  proposed,  based  on  the  Early- 
Late  method.  This  bit  synchronization  scheme  then  led  to  the  postulation 
of  a  complete  receiver  algortihm  including  interference-tracking, 
phase  estimation,  and  bit  synchronization.  A  block  diagram  of  the 
algorithm  was  given. 

Task  8. 

The  Monte  Carlo  simulation  routine  used  and  reported  previously 
[1,  4]  was  restructured  and  re-written.  The  routine  was  simplified  con- 
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siderably  and  was  augmented  to  accomodate  the  new  detection  and  phase¬ 
tracking  algorithms.  The  chief  reason  for  this  effort  was  to  achieve 
shorter  run  times  in  line  with  restrictions  imposed  by  the  ASO  Computer 
Facility  (CDC-6600). 

The  performance  of  the  IDEI  detector  with  phase  tracking  was  eval¬ 
uated.  It  was  found  that  the  performance  was  considerabley  degraded  over 
previous  results  for  perfect  phase  references.  Two  possible  causes  for 
the  degradation  were  discussed. 

In  summary,  further  research  on  the  IDEI  algorithms  is  recommended 
in  the  following  areas.  Most  importantly,  a  method  of  phase  estimation 
should  be  sought  wherein  the  phase  estimator  is  imbedded  in  the  inter¬ 
ference  tracking  filter.  The  purpose  is  to  reduce  the  effects  of  the 
large  additive  colored  noise  upon  the  phase  estimator.  Rather  than 
tracking  phase  in  parallel  with  the  colored  noise  tracking  filters,  phase 
should  be  tracked  after  the  colored  noise  has  been  removed  from  the 
data.  Secondly,  further  effort  should  be  devoted  to  optimizing  the 
multiplicative  nois\s  tracking  filter  for  the  non-Gaussian  perturbations 
produced  by  the  phase  variations.  Finally,  the  proposed  bit  synchroni¬ 
zation  algorithm  should  be  studied  and  evaluated. 


APPENDIX  A 


THE  CLOSED-FORM  ERROR- 
RATE  PROGRAM 


(This  appendix  contains  listings  of  the  newly  written  simulation 
program  and  the  closed-form  error-rate  evaluation  program 
reported  previously.) 


C  THIS  IS  MAIN  PROGRAM  FOR  THE  CLOSED-FORM  ERROR  RATE  FOR 
C  IMPERFECT  IDENTIFICATION  WHICH  IS  A  EXTENSION  OF  PROGRAM  YOONM5 
C  **  REQUIRED  SUBROUTINE  ** 

C  <1>  RES3  i  MAIN, DATA,  INPUT1,  INPUT2, PARALL,  PREPAR 
C  (2)  CFERAT  ;  CFERAT, WKFLT, ERF 
C  (3)  VTT  i  VTT.  CAYLEY, QAUS 

C  < 4 )  EIGEN 
C  <  5 )  COMAT 
C  REMARK 

C  (1)  CHECKING  THE  CLOSED-FORM  ERROR  RATE  FOR  PERFECT 

C  IDENTIFICATION,  SET  IMODE  1  AVOIDING  THE  SAME  EIGEN-VALUE 

C  IN  SUB.  CAYLEY 

C  (2)  TO  GET  THE  STEADY-STATE  KALMAN  GAIN.  SET  KSMAX  50-100 

C  IN  GENERAL 

C  (3)  ESTIMATED  TRANSITION  MATRIX  PHEER  AND  DPHEE  ARE  VARYED 

C  IN  SUBROUTINE  INPUT1  AND  ESTIMATED  KALMAN  GAIN  GSTAR  IS 

C  VARYED  IN  SUBROUTINE  INPUT2  EACH  TIME. 

C  PROGRAMMER 
C  CHANG-JUNE  YOON 

C  ELECTRICAL  ENGINEERING  DEPT. 

C  TEXAS  A  &  M  UNIVERSITY 

C 

COMMON/ORDER/N.  N2 
COMMON/SAMPLE/NSPB, TB, TBR 
COMMON/OPTION/NOS,  AEST 

COMMON/R AT I O/ENODB , ENODBR. SJRDB, SJRDBR 
COMMON/GDB/GN.  GNR,  GJ,  GJR 
COMMON /WORNOW/ IMODE,  KSMAX,  IOJ 
COMMON/FREQ/FZ.  FP ( 3 ) 

COMMON/PAR AM/OAMMA ( 6,  2 ) , PHEE < 6, 6 ) , H< 2, 6), Q(2.  2), R<2,  2) 
COMMON/PARAMR/PHEER < 6,  6), DPHEE <6, 6) , GSTAR <6, 2 ) ,  BSTAR <2,  2) 

CALL  ASSIGN<5,  'SY:  RES3.  DAT ',  1 1 ,  'RDO'»  'NC',1> 

C 

CALL  DATA 

C 

C  NOP TNI 

C  1 , NO  CHANGE 

C  2, CHANGE  ENODB 

C  3, CHANGE  SJRDBR 

C  4, CHANGE  NSPB, GK 

C  N0PTN2 

C  1 . NO  CHANGE 

C  2, CHANGE  ENODB 

C  3, CHANGE  SJRDBR 

C  4. CHANGE  GK 

C  5, CHANGE  SJRDB, SJRDBR 

C  IF  NOP TNI, N0PTN2  IS  1,  THEN  NCASE1, NCASE2  IS  1  RESPECTIVELY 

C  NCASE1  ;  NUMBER  OF  CASE  FOR  N0PTN1 

C  NCASE2  j  NUMBER  OF  CASE  FOR  N0PTN2 

C  I PAR AM 

C  0. NO  PRINT-OUT  PARAMETERS  AND  STATISTICS  IN  INPUT1  AND  INPUT2 

C  1, PRINT-OUT 

C  IOV 

C  0,  NO  CALCULATION  KALMAN  GAIN  FOR  A  CORRECT  PARAMETERS  INPUT1. 

C  1, CALCULATION. 

C 

READ <  5, 701 >  N0PTN1, N0PTN2.  NCASE1,  NCASE2,  IPARAM,  IGV 

701  FORMAT  <  61 5 ) 

READ (5, 702)  GK 

702  FORMAT (El 5.  6) 

DO  2000  11=1, NCASE1 
00  TO  (1, 2,  3,  4),N0PTN1 

1  00  TO  50 

2  READ (5, 70S)  ENODB 
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ENODBR * t NUDB 
GO  TO  90 

3  READ (9,  705)  SJRDBR 
GO  TO  90 

4  READ <9,  707)  N8PB,  GK 
90  CONTINUE 

709  FORMAT (El 5.  6) 

706  FORMAT (2E 15.  6) 

707  FORMAT  (1 5,  El  9.  6) 

C 

DO  1000  111=1,  NCASE2 

GO  TO  (11.  12,  13,  14,  15>,N0PTN2 

11  GO  TO  60 

12  READ (5, 705)  ENODB 
ENODBR-ENODB 

00  TO  60 

13  READ(5. 705)  SJRDBR 
GO  TO  60 

14  READ (9. 705)  GK 
GO  TO  60 

15  READ(5, 706)  S JR DB, SJRDBR 
60  CONTINUE 

C 

WRITE (6.  650)  NSPB,  TB,  ENODB,  SJRDB,  SJRDBR,  GK, AEST 
650  FORMAT (2X»  5HNSPB=,  15,  2X,  3HTB*.  E13.  6,  2X,  6HEN0DB=,  E13.  6,  2X, 

16HSJRDB-,  E13.  6.  2X.  7HSJRDBR=»  E13.  6,  2X,  3HGK=,  E13.  6,  2X,  5HAEST=»  E13.  6) 
CALL  INPUTKIPARAM.  IGV) 

CALL  INPUT2( IPARAM. OK) 

C 

CALL  CFER AT ( ERATCL ) 

C 

WRITE (6, 651)  ERATCL 

691  FORMAT (2X.26HCL0SED-F0RM  ERROR  RATE  IS  , 30X,  10H*********=.  E13.  6) 
1000  CONTINUE 

WRITE (6,  751) 

791  FORMAT! 5X.  11HEND  OF  CASE,//) 

2000  CONTINUE 
STOP 
END 

C 

SUBROUTINE  DATA 

COMMON/ORDER /N,  N2 

COMMON /SAMPLE /NSPB,  TB,  TBR 

COMMON/OPTION/NOS.  AEST 

COMMON/RATIO/ENODB,  ENODBR.  SJRDB, SJRDBR 

COMMON/WORNOW/IMODE,  KSMAX,  IOJ 

COMMON/FREQ/FZ.  FP(3) 

C  N, N2  :  SYSTEM  ORDER 

C  NOS  :  (1)  PSK,  (2)  F8K 

C  AEST  :  SIGNAL  MAGNITUDE  IN  SUB.  REFOEN 

C  IMODE:  (1)  DIAGONAL  PHEE  MATRIX  AND  PERFECT  IDENTIFICATION 
C  (2)  DIAGONAL  PHEE  MATRIX 

C  (3)  GENERAL  IMPERFECT  IDENTIFICATION 

C  KSMAX:  MAXIMUM  NUMBER  OF  ITERATION  FOR  8TE ADY-ST ATE  KALMAN  CAIN 
C  IGV  :  (0)  NO  CALCULATION  FOR  CORRECT  KALMAN  GAIN  AND  VINOV  IN  INPUT1 
C  (1)  CALCULATION  FOR  CORRECT  KALMAN  GAIN  AND  VINOV  IN  INPUT1 

C  FZ.FP:  ZERO, POLE  FREQUENCY  FOR  LOW-PASS  FILTER 
READ!  9,  600)  N.  N2 
READ (5. 601)  NSPB, TB. TBR 
READ(5,  602)  NOS, AEST 

READ ( 9. 603)  ENODB, ENODBR , SJRDB, SJRDBR 
READ(5. 604)  IMODE, KSMAX, IOJ 
READ(9,  603)  FZ,  (FP(  I ),  1  =  1, 3) 

600  FORMAT (21 9) 

601  FORMAT  (15,  2E19.  6) 

602  FORMAT ( 19,  El 9.  6) 
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“'503  FCBfHAT  f4ET571D  ~  - -  - - 

604  FORMAT (31 5) 

RETURN 

END 

SUBROUTINE  INPUT1 ( IPARAM, IGV) 

TO  QET  THE  REAL  PARAMETERS  AND  STATISTICS  GIVEN  VALUES. 

ALL  WE  NEED  IN  HERE  ARE  GAMMA.  PHEE.  Hi  R 
GAIN  AND  VI NOV  ARE  FOR  REFERENCE 

IF  IGV  :  0  -  NO  CALCULATION  FOR  CORRECT  KALMAN  GAIN  AND  VINOV 
:  1  -  CALCULATION  FOR  CORRECT  KALMAN  GAIN  AND  VINOV 

THEREFORE  OK  ALWAYS  SET  1.  .  FOR  PERFECT  IDENTIFICATION. 

COMMON /OR DER/N.  N2 
COMMON/SAMPLE/NSPB. TB. TBR 
COMMON/OPTION/NOS, AEST 

COMMON/RAT IO/ENODB , ENODBR. SJRDB. SJRDBR 
COMMON/GDB/GN. GNR. GJ. GJR 
COMMON/UORNQW/ IMODE.  KSMAX.  IOJ 
COMMON/FREG/FZ.  FP(3> 

COMMON /PAR AM /GAMMA ( 6. 2).  PHEE<6, 6). H(2, 6), Q<2. 2), R<2. 2) 
DIMENSION  VINQV<2.  2).  GAIN(6. 2) 

CALL  PARALL(  1.  .  BN. GAMMA. PHEE, H, ENODB. SJRDB, GN, GJ, R,  IGV 
1. GAIN, VINOV) 

IF( IPARAM.  EG. 0)  GO  TO  40 
WRITE (6,  610) 

610  FORMAT  <2X,  32H*REAL  PARAMETERS  AND  STATISTICS*,/) 

DO  20  1=1, N 

WRITE<6, 611)  I ,  GAMMA ( 1 ,  1),  I ,  PHEE  < I »  I),  I,  H(l,  I) 

611  FORMAT  (2X.  5HGAMD  < .  II,  2H)=,  E13.  6,  2X,  5HPHIDC,  II,  2H)=,  E13.  6, 

12X,  3HHT ( ,  II.  2H)=,  E13.  6) 

20  CONTINUE 

IF< IGV.  EG. O)  GO  TO  40 
WRITE(6, 615)  GN 
615  FORMAT < /, 2X.  3HGN=, E13.  6) 

WRITE<  6, 612) 

612  FORMAT </.  2X,  5HGAIN=,  26X.  6HVIN0V*) 

DO  25  1=1, N2 

IF <  I .  GT.  2)  GO  TO  30 

WRITE <6,  613)  (GAINd,  J),  J=l,  2),  (VINOVO,  J).  J=l,2> 

613  FORMAT ( 2X,  2E 13.  6.  5X.  2E13.  6) 

GO  TO  25 

30  WRITE<6,  614)  (GAINd,  J),  J=l,  2) 

614  FORMAT  ( 2X ,  2E 1 3.  6 ) 

25  CONTINUE 

40  CONTINUE 

WRITE(6, 617)  BN 

617  FORMAT (/,  2X,  2 1 HEQU I VALENT  BANDWIDTH*, E13. 6, /> 

RETURN 

END 

SUBROUTINE  INPUT2 ( IPARAM, GK ) 

THIS  SUBROUTINE  GSTAR  AND  DPHEE  FOR  DIFFERENT  FILTER  BANDWIDTH 
THESE  GSTAR  AND  DPHEE  WITH  GAMMA,  PHEE, H, R  ARE  USED  TO  CALCULATE 
RESIDUAL  VARIANCE  IN  SUBROUTINE  CFERAT  AND  VTT. 

THEREFORE  IGV  ALWAYS  SET  1  HERE. 

COMMON/ORDER /N,  N2 

COMMON/RATIO/ENODB, ENODBR, SJRDB, SJRDBR 
COMMON/GDB/GN,  ONR, GJ, GJR 
COMMON/WORNOW/ IMODE,  KSMAX,  IOJ 

COMMON/PAR AM/GAMMA ( 6,  2), PHEE(6, 6), H(2, 6), G(2, 2), R(2, 2) 
COMMON/PAR AMR /PHEER (6,  6),  DPHEE (6,  6),  GSTAR (6, 2) , BSTAR (2, 2) 
DIMENSION  GAMMAR ( 6, 2), RR<2, 2), VIN0VR(2. 2) 

CALL  PARALL(GK, BNR, GAMMAR, PHEER, H, ENODBR, SJRDBR, GNR, OJR,  RR 
1.1,  GSTAR,  VINOVR) 

DO  10  1=1. N2 
DO  10  J-1.N2 


10  DPHEE (I .  J ) “PHEE (I ,  J ) -PHEER  (I .  J ) 

IF  (IP  ARAM.  EG.  0)  RETURN 
WRITE<6.  600) 

600  FORMAT </,  2X,  39HEST I MATED  PARAMETERS  AND  STATISTICS* / ) 

DO  20  I-l.N 

WRITE<6.  601)  I  ,  GAMMAR  ( I  .  1),  I.  PHEER  (I,  I),  I,H<1,  I),  I,  DPHEE  <  I,  I) 

1,  I,  OSTAR<  I,  1 ) 

601  FORMAT  ( 2X  *  '  OAMDR  <  '*  II*  ')«'*E13.  6,  2X,  'PHIDR  (  11*  ')=',E13.  6 

1, 2X,  'HTR  (  II.  ')='.  E13.  6,  5X.  'DPHEE (  II,  ')-',  E13.  6,  2X,  '©STAR <  ' 

2,  II.  ')=',  E13.  6) 

20  CONTINUE 

WRITE<6,  602)  BNR 

602  FORMAT ( /,  2X, 21HEGUI VALENT  BANDWIDTH-, E13.  6) 

WRITEI6. 603)  VINOVR(l.l) 

603  FORMAT < / .  2X*  12HVIN0VR(1*  1)-, E13.  6) 

RETURN 

END 

C 

SUBROUTINE  PARALL(OK* BN.  GAMMA*  PHEE, H*  ENODB. SJRDB,  GN. GJ, R,  IGV 
1.  GAIN.  VINOV) 

COMMON/ORDER/N,  N2 
COMMON/OPTION/NOS,  AEST 
COMMON/SAMPLE/NSPB.  TB*  TBR 
COMMON/ WORNOW/ 1 MODE* KSMAX,  IOJ 
COMMON/FREQ/FZ* FP(3) 

DIMENSION  FPR ( 3 ) 

DIMENSION  GAMD(3),  PHID(3).  HT(3) 

DIMENSION  GAMMA ( 6, 2),  PHEE<6,  6) . H(2* 6). Q(2, 2), R<2, 2) 

DIMENSION  PVPT  <  6.  6),  GTG(6.  6) ,  VEST (6,  6) ,  VPRED< 6,  6 ) .  HVHT<2, 2) 

1,  VINOV< 2,  2),  VINV<2.  2),  VPHT(6,  2).  GAIN(6* 2), GH(6.  6) 

REAL  IMGH(6,  6) 

PI-4.  *ATAN(  1.  ) 

DELPHI-.  785 

DELMEO-DELPH I *2 .  *P I /TB 
IF  (NOS.  EQ.  1)  GO  TO  1 
SUMF-O. O 
DO  2  K-l.NSPB 

2  SUMF— SUMF+(SIN( (K— .  5)*TB*DELMEG/NSPB) )**2 
CONSTF -SORT ( SUMF ) 

ON— C0NSTF*10.  ** (— EN0DB/20.  ) 

GO  TO  3 

1  CONSTP— SORT (NSPB/2.  ) *ABS( SIN ( DELPHI ) ) 

ON— C0N8TP*10.  **  <  — EN0DB/20.  ) 

3  0J-10.  **(-SJRDB/20.  )/S0RT<2.  ) 

R ( 1 . 1 ) —ON** 2 

R  ( 1 , 2 )  — O.  O 
R<2.  1  )-0.  0 
R (2. 2  >  — 0N**2 
FZR— OK*FZ 
DO  9  I-l.N 
9  FPR ( I ) — QK*FP ( I ) 

T— T1/N8PB 

CALL  PREPAR(T,  FZR*  FPR.  QAMD.  PHID,  HT, BN.  IOJ) 

DO  10  I-1.N2 
DO  10  J— 1,  2 

10  GAMMA  < I ,  J ) -0.  0 
DO  11  I-l.N 
GAMMA(I.  1 ) — GAMD (I ) 

11  GAMMA < I+N.  2 ) -GAMD ( I ) 

C  NEW  WEIGHTED  GAMMA  MATRIX 
DO  12  I-1.N2 
DO  12  J-1,2 

12  OAMMAd.  J)-CJ*OAMMAd»  J) 

DO  19  I-1.N2 

DO  19  J-l, N2 
15  PHEEd,  J)-0.  0 
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DO  ~\6  I=T7R  - -  - -  — 

PHEEd, I )=PHID( I ) 

16  PHEEd+N,  I+N)=PHID!I) 

DO  20  1=1,  2 
DO  20  J=1.N2 

20  H( I. J>=0. O 
DO  21  1=1. N 
H< 1 »  I ) =HT  < I ) 

21  H!2,  I+N)=HT < I ) 

C 

IF (  IGV.  EQ.  O)  RETURN 

C  CALCULATE  THE  STEADY-STATE  KALMAN-GAIN 
DO  32  1  =  1.  N2 
DO  32  U=1.N2 
32  VEST < I » J)=0.  O 
DO  35  KS=1 »  KSMAX 

CALL  MABCT <PHEE.  N2,  N2.  VEST.  N2.  PHEE.  N2.  PVPT,  6,  6.  6.  6.  6.  6,6,6) 

CALL  MATMUL  !  2.  GAMMA,  N2,  2,  GAMMA,  N2,  GTG.  6,  2,  6.  2,  6.  6) 

CALL  MATAS <1, PVPT, N2, N2, GTG, VPRED, 6, 6) 

CALL  MABCT  (H,  2,  N2,  VPRED,  N2,  H,  2,  HVHT,  2.  6,  6,  6,  2,  6,  2,  2) 

CALL  MATAS(1, R,  2,  2, HVHT,  VINOV,  2,  2) 

CALL  MATMUL  (2,  VPRED,  N2,  N2,  H,  2,  VPHT,  6,  6,  2,  6,  6,  2) 

DET =V  I  NOV  <  1 »  1)*VIN0V!2,  2>— VINOVd,  2>*VIN0V<2,  1) 

VINV(1, 1 )=VIN0V<2» 2) /DET 
VINVC 1, 2)=— VINOV< 1,2) /DET 
VINV(2,  1 )=— VIN0V(2,  1>/DET 
VINV<2, 2 ) = V I NOV  < 1 .  1)/DET 

CALL  MATMUL  d,  VPHT,  N2,  2,  VINV,  2,  GAIN,  6,  2,  2,  2,  6,  2) 

CALL  MATMUL (1,  GAIN,  N2,  2,  H,  N2,  GH,  6,  2,  2,  6,  6,  6) 

DO  36  1=1.  N2 
DO  36  J=1,N2 
IMGHd,  J )  =— GH ( I ,  J) 

IFd.EG.  J)  IMGHd,  J)  =  l.  0-GH(I,  J) 

36  CONTINUE 

CALL  MATMUL  <  1,  IMGH,  N2,  N2,  VPRED,  N2,  VEST,  6,  6,  6.  6,  6,  6) 

35  CONTINUE 
RETURN 
END 

SUBROUTINE  PREPARd,  FZ,  FP,  GAMD,  PHID,  HT,  BN,  INOPT) 

C  *-»**•»■»-*****•»■»********■«-********##*****##**###***#*********#*#****## 

C  PREPAR:  MODIFICATION  SUBROUTINE  ADDED  TO  SUBROUTINE  INPUT 

C  TO  PERFORM  PRE-CALCULATIONS  OF  FILTER  PARAMETERS 

C  *1/0  PARAMETERS* 

C  *  INPUT  * 

C  T:  SAMPLING  TIME 

C  FZ:  ZERO  FREQUENCY 

C  FP:  POLE  PREQUENCIES  <3) 

C  INOPT:  1-DIGIT  CODE  FOR  SELECTION  OF  REAL/COMPLEX  ZERO 

C  AND  UNITY  GAIN/VARIANCE  FOR  FILTER  PARAMETER  CALCULATIONS 

C  =1, REAL  ZERO,  UNIT  GAIN 

C  =2, REAL  ZERO,  UNIT  VARIANCE 

C  =3, COMPLEX  ZERO, UNIT  GAIN 

C  -4, COMPLEX  ZERO, UNIT  VARIANCE 

C  *  OUTPUT  * 

C  PHID:  FILTER  TRANSITION  WEI0HTS!3> 

C  GAMD:  FILTER  INPUT  WEIGHTS! 3) 

C  HT:  FILTER  OUTPUT  WEI0HTS<3> 

C  BN:  EQUIVALENCE  NOISE  BANDWIDTH 

C  *  INTERNAL  FILTER  PARAMETERS  * 

C  Z:  ZERO  FREQUENCY,  IN  RADIANS 

C  P:  POLE  FREQUENCIES  (3).  IN  RADIANS 

C  R:  RESIDUES! 3) 

C  RE:  RESIDUES!3) 

C  OAINK  GAIN  CONSTANT 

C  #*******************************************#*##*****#**' 

DIMENSION  FP!3),P!3)»R!3),  RE! 3) ,  PHID! 3) ,  0AMD!3) ,  HT!3) 
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n  n  o 


PI-4.  *AIAN<  1.  ) 

IF ( INOPT.  GT.  2)  00  TO  100 
C  FREQUENCY  CALCULATIONS 
Z»<-2.  )*PI#FZ 
DO  1  1-1.3 

1  P< I >  —  <-2.  >*PI*FP ( I ) 

C  RESIDUE  CALCULATIONS 

DO  5  1-1.3 
D— 1. 

DO  10  J-l.  3 

IF<  I.  EQ.  J)  00  TO  10 

D— D*  <P<D— P<J) ) 

10  CONTINUE 

R(I)« <P(I )— Z)/D 
5  CONTINUE 

C  TRANSITION  WEIGHTS 

DO  20  1-1.3 
20  PHID( I )— EXP(P( I  )*T) 

C  INPUT  WEIGHTS 

DO  29  1-1.3 

25  0AMD(I)*(1.  — PHID( I ) )*R< I >/<— P< I ) ) 

C  UNITY  GAIN 

IF<  INOPT.  NE.  1)  GO  TO  30 
OAINK— P< 1 )*P(2)*P(3)/Z 
GO  TO  35 
30  CONTINUE 
C  UNITY  VARIANCE 

SUM-0.  0 
DO  40  1-1.3 
DO  40  J-l.  3 

40  SUM— SUM+OAMD ( I ) *C  AMD <  J ) / ( 1 .  0-PHID< I >*PHID( J) ) 

OAINK— 1 .  /SORT (SUM) 

35  CONTINUE 

C  NEW  WEIGHTED  INPUT  MATRIX 
DO  36  1-1.3 

36  OAMD ( I ) — OA INK*CAMD ( I ) 

C  OUTPUT  WEIGHTS 

DO  45  1-1.3 
45  HT  < I )  —  1 . 

C  EQUIVALENT  NOISE  BANDWIDTH 

DO  50  1-1.3 
D-l. 

DO  55  J-l,  3 

IF< I .  EQ.  J)  00  TO  55 

D-D* ( P  < I ) **2-P  <  J ) **2 ) 

55  CONTINUE 

RE (I )-0AINK**2*<P( I >**2-Z**2)/<2.  *P( I )*D) 

50  CONTINUE 

00-GAINK*Z/(P(l)*P(2)*P(3) ) 

BN—  < RE ( 1 ) +RE ( 2 ) +RE ( 3 ) ) / ( 2.  *00**2) 

RETURN 
100  CONTINUE 

MODIFIED  TRAN8FER  FUNCTION  HAVING  COMPLEX  ZERO. 

FREQUENCY  CALCULATION. 

Z**2-P < 2 ) **2-2*P < 1 ) **2>  TO  HAVE  A  JW-AXIS  ZERO, Z  SHOULD  BE  POSITIVE 
Z— 2.  *PI*FZ 
P<l)-Z 

P<2)-8QRT (3.  )*P(1) 

P<3)— 2.  *PI*FP<3> 

FP <  1 )  — P ( 1  > / < -2.  *PI> 

FP ( 2 )  — P < 2 ) / <  —2.  *PI> 

FP< 3)-P<3)/<-2.  *PI > 

C  RESIDUE  CALCULATIONS 
DO  110  1-1,3 

D-l.  46 

DO  120  J-l.  3 


gf>'  r  I  kai'if 


Iff  I.  EG  JTOtr  TD  120 - 

D— D*  <  P  < I ) — P ( J ) ) 

120  CONTINUE 

R<I)-<P<I)**2+Z**2)/D 
110  CONTINUE 

C  TRANSITION  WEIGHTS 
00  125  1  =  1,3 
125  PHID< I )=EXP<P( I  )*T) 

C  INPUT  WEIGHTS 

DO  130  1  =  1,  3 

130  GAMD ( I )  —  <  1 .  — PHID( I > )#R( I )/( 1.  O-PHIDC I )*PHID( J) ) 

C  UNITY  GAIN 

IF( INOPT.  NE.  3)  GO  TO  135 
GAINK=2.  *P I *FP  < 1 >  *FP  <  2  >  *FP ( 3 ) /FZ**2 
00  TO  140 
135  CONTINUE 
C  UNITY  VARIANCE 

SUM=0.  0 
DO  150  1=1, 3 
DO  150  J-1.3 

1 50  SUM-SUM+CAMD ( I ) *GAMD ( J) / ( 1 .  O-PHID ( I ) »PHID(J) ) 
GAINK=1.  /SORT (SUM) 

140  CONTINUE 

C  NEW  WEIGHTED  INPUT  MARTRIX 
DO  141  1=1,3 

141  GAMD ( I ) =GA I NK*QAMD ( I ) 

C  OUTPUT  WEIGHT 

DO  155  1  =  1,3 
155  HT ( I )  =  1 . 

C  EG I VALENT  NOISE  BANDWIDTH 

DO  160  1=1,  3 
D=1 . 

DO  170  U-1 , 3 

IF ( I .  EQ. J)  GO  TO  170 

D-D*  <  P ( I ) **2-P  <  J ) **2 ) 

170  CONTINUE 

RE( I )  —  <  — 1 .  )*GAINK**2*(P(I)**2+Z**2)**2/(2.  *P(I)*D) 
160  CONTINUE 

G0= ( — 1 . )*GAINK*Z**2/<P< 1 )*P(2)*PC3) ) 

CN— ( RE ( 1 ) +RE < 2 ) +RE < 3 ) )/<2.  *00**2) 

RETURN 

END 
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SUBROUTINE  CFERAT ( ERATCL ) 

EXTERNAL  ERF 
COMMON/ORDER/N, N2 
COMMON /SAMPLE /NSPB, TB, TBR 
COMMON/ WORNOW/IMODE,  KSMAX,  IOJ 

COMMON/PARAM/GAMMA ( 6,  2),  PHEE<6,  6),  H(2,  6),  Q<2»  2)>  R!2,  2) 
COMMON/PARAMR/PHEER ( &, 6) ,  DPHEE < 6/  6),  GSTAR! 6,  2),  BSTAR(2, 2) 
DIMENSION  XEST1 (6) » XEST2(6), XPRED1 <6> » XPRED2<6) 

DIMENSION  SIG1!2>,  SIG2!2),  ESI (2).  ES2(2) 

DIMENSION  B1 (2, 300), VTTJ!2, 2),  VTILDA!300) 

DIMENSION  VXX  <6,  6), TEMPI6,  6).  GTG!6,  6),  F<6,  6).  VXXT!A,  6), VXXT1I6.  6) 

1,  VXXT2!6,  6 ) »  VXTXT  ( 6,  6).  VXTXT1  <6,  6).  VXTXT2I6.  6),  VXTXT3!6.  6) 

2,  VXTXT4!6,  6), VXTXT5(6,  6),  GRG!6,  6),  TEMPI (6.  6) 

REAL  IMGHI6, 6) 

CALL  MATMUL  ( 2,  GAMMA,  N2,  2,  GAMMA,  N2,  GTG,  6,  2.  6.  2,6,6) 

CALL  MATMUL ( 1,  GSTAR,  N2.  2,  H,  N2,  TEMP,  6,  2,  2,  6,6,6) 

DO  5  1  =  1.  N2 
DO  5  J=1.N2 
IMCH(  I ,  J)— TEMPI  I,  J) 

IF!  I.  EG.  J)  IMGH<  I,  J)  =  l.  O-TEMP  <  I,  J) 

5  CONTINUE 

CALL  MATMUL!  1,  PHEER,  N2,  N2,  IMGH,  N2,  F,  6.  6,  6,  6,6,6) 

INITIALIZE  VXX, VXXT  AND  VXTXT 
DO  6  1  =  1,  N2 
DO  6  J= 1,  N2 
VXX(I.  J)=0.  0 
VXXT <  I,  J)=0.  O 

6  VXTXT ( I, J)=0. 0 

DO  1  KS=1, KSMAX 

VXX (K)  =  PHEE  *  VXX(K-l)  *  PHEE '  +  GAMMA  *  G  *  GAMMA' 

CALL  MABCTIPHEE.  N2,  N2,  VXX,  N2,  PHEE,  N2,  TEMP.  6.  6,  6,  6.  6,  6,6.6) 

CALL  MATAS! 1,  TEMP.  N2.  N2,  GTG,  VXX,  6,  6) 

VXXT !K ! K— 1 )  =  PHEE  *  VXXT<K-1 ! K-2)  *  F'  +  PHEE  *  VXXIK)  *  DPHEE' 

+  GAMMA  *  Q  *  GAMMA' 

CALL  MABCTIPHEE.  N2,  N2,  VXXT,  N2.  F,  N2,  VXXT1,  6,  6,  6,6,  t,  6,  6,  6) 

CALL  MABCT! PHEE,  N2,  N2,  VXX,  N2,  DPHEE,  N2,  VXXT2,  6,  6,  6.  6,  6,  6,  6,  6) 

CALL  MATAS! 1.  VXXT1,  N2,  N2,  VXXT2,  TEMP,  6, 6) 

CALL  MATAS! 1, TEMP.  N2,  N2, GTG, VXXT, 6,  6) 

VXTXT !K+1 ! K)  =  F  *  VXTXT!K!K-1)  *  F'  +  2.  *  DPHEE  *  VXXT!KfK-l)  *  F' 
+  DPHEE  *  VXX ! K )  *  DPHEE'  +  GAMMA  *  Q  *  GAMMA' 

+  PHEER  *  GSTAR  *  R  *  GSTAR'  *  PHEER ' 

CALL  MABCT  IF,  N2,  N2,  VXTXT,  N2,  F,  N2,  VXTXT1,  6,  6,  6,  6,  6,  6,  6,  6) 

CALL  MABCT! DPHEE,  N2.  N2,  VXXT,  N2,  F,  N2,  VXTXT2,  6,  6,  6,  6,  6,  6,6,6) 

DO  IS  1=1,  N2 
DO  13  J=1,N2 

13  VXTXT3! I,  J)=VXTXT2! J,  I ) 

CALL  MABCT! DPHEE,  N2,  N2,  VXX,  N2,  DPHEE,  N2,  VXTXT4,  6,  6,  6,  6,  6,  6,  6,  6) 
CALL  MABCT! GSTAR,  N2,  2,  R,  2,  GSTAR,  N2,  ORG.  6,  2,  2,  2,  6,  2,  6,  6) 

CALL  MABCT ! PHEER,  N2,  N2,  GRG,  N2,  PHEER,  N2,  VXTXT5,  6,  6,  6,  6.  6,  6,  6.  6) 
CALL  MATAS! 1, VXTXT1,  N2,  N2, VXTXT2,  TEMP, 6, 6) 

CALL  MATAS! 1.  TEMP.  N2,  N2,  VXTXT3,  TEMPI, 6, 6) 

CALL  MATAS! 1,  TEMPI,  N2,  N2,  VXTXT4, TEMP, 6, 6) 

CALL  MATAS! 1,  TEMP,  N2,  N2,  GTG,  TEMPI,  6,  6) 

CALL  MATAS! 1,  TEMPI,  N2,  N2,  VXTXT5,  VXTXT, 6,  6) 

1  CONTINUE 

DO  25  1  =  1,  N2 
XEST1 ! I )=0.  0 
XE8T2! I >=0.  0 
25  CONTINUE 
C  SIG1  ES!M=0, N=0) 


C  5102  ESTH=?57R»1') - - -  - 

A=0  O 

DO  20  K-l.NSPB 

CALL  REFOEN (K, O, FTRO, GTRO,  FRRO, ORRO) 

CALL  REFOEN ( K,  1 . FTR 1 , OTR 1 . FRR i ,  ORR 1 ) 

SI01 < 1 )"FTRO-FRRO 
SI01 <  2 ) -OTRO-ORRO 
SIG2( 1 ) “FTRO— FRR 1 
S I 02 ( 2 ) —GTRO-GRR 1 

CALL  WKFLT(K,  XEST1,  XPRED1 ,  SI01 , ESI ) 

CALL  WKFLT <K# XEST2. XPRED2. SI02, ES2 ) 

A-A+  <  ES2  < 1 ) **2+ES2 ( 2 ) **2-ES 1 < 1 > **2-ES 1(2) **2  > 

B1  ( 1,  K)**ES1  <  1  )-ES2(  1 ) 

B 1 <  2, K ) =ES 1(2) -ES2 ( 2 ) 

20  CONTINUE 

C  THE  VII.  INNOVATION  VARIANCE,  IS  DECOUPLED,  SO  IS  VTTJ  SINCE 
C  THE  TEST  SYSTEM  IS  DECOUPLED. 

C  IF  THE  SYSTEM  IS  COUPLED,  THEN  THE  EVALUATION  OF  MEAN  AND  VARIANCE 
C  MUST  BE  MODIFIED. 

DO  30  J=1 , NSPB 
L=J-1 

CALL  VTT(L, F, VXTXT,  VXXT,  VTTJ) 

VT I LDA ( J ) =VTT J ( 1 , 1 ) 

30  CONTINUE 
B-0.  O 

DO  35  J*1 , NSPB 
DO  35  K=1 « NSPB 
L= I ABS ( J— K ) + 1 

B=B+(B1 (1, J ) *B 1 ( 1 , K ) +B 1(2,  J)*B1 <2,  K) )*VTILDA(L> 

35  CONTINUE 

IF ( B.  LE.  O.  )  GO  TO  50 
SIGMAB=SQRT(B) 

X=A/SIGMAB 

SUFX—X/ (2.  *SQRT(2.  )) 

ERATCL=0.  5*( 1.  -ERF(SUFX)) 

WRITE (6, 600)  A, SIGMAB, VTILDA( 1 )» X 

600  FORMAT (2X,  'MU=',E13.  6,  2X,  'SICMA=  ' ,  E13.  6,  2X,  'VTT(0)= ',  E13.  6. 

12X.  /MU/SIGMA=',  E13.  6) 

RETURN 

50  WRITE (6. 601 ) 

601  FORMAT ( 2X, 20HVAR I ANCE  IS  NEGATIVE) 

DO  60  I=»l ,  NSPB 

60  WRITE (6, 602)  I , VTILDA( I ) , B 1 < 1 , I ) , B 1 ( 2, I ) 

602  FORMAT (2X» 7HVTILDA(,  13, 2H)=, E13.  6, 2X,  15HTRACKING  ERROR-,  2E15.  6) 
RE1URN 

END 

SUBROUTINE  WKFLT (KS,  XEST,  XPRED,  SIG,  V) 

COMMON/ORDER /N,  N2 

COMMON/PARAM/GAMMA ( 6,  2), PHEE(6,  6) , H(2,  6 ) , Q<2, 2),  R(2,  2) 
COMMON/PAR AMR/PHEER ( 6,  6), DPHEE(6, 6 ) ,  GSTAR ( 6,  2 ) ,  BSTAR (2,  2) 
DIMENSION  XEST (6) , XPRED(6), SI0(2), V(2), ZHAT(2>,  GV(6) 

CALL  MATVEC ( PHEER , N2,  N2, XEST,  XPRED,  6,  6) 

CALL  MATVEC (H.  2,  N2, XPRED,  ZHAT,  2,  6) 

CALL  VECAS(2, SIG, ZHAT, V,  2) 

CALL  MATVEC (GSTAR, N2,  2,  V,  GV, 6,  2) 

CALL  VECAS(1, XPRED, OV,  XEST,  6) 

RETURN 

END 

SUBROUTINE  REFOEN (KS. M, FTR, OTR, FRR, ORR) 

COMMON/SAMPLE/NSPB,  TB,  TBR 
COMMON/OPTION/NOS,  AEST 
TK—(KS— 0. 5) /NSPB 
TKRMOD- ( TK- I F I X ( TK ) ) *TBR 
DELPHI*  785 

DELME0-DELPHI*8.  *ATAN(1.  )/TB 
IF  (NOS.  NE.  1)  00  TO  l 


49 
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IF<n.  EQ.  0)  PHfct IH-DfcLPHl 
IF<M.  EQ.  1 )  PHEETR— DELPHI 
00  TO  2 

1  IF<M.  EQ.  0)  PHEETR«DELMEQ*TKRMOD 
IF<M.  EQ.  1>  PHEETR— DELMEO*TKRMOD 

2  FTR«COS<P*«ETR> 

©TR-SIN< PHEETR) 

FRR-AEST*COS ( PHEETR  > 

©RR“AEST*SIN<PHEETR ) 

RETURN 

END 

FUNCTION  ERF(X) 

THIS  IS  AN  APPROXIMATION  OF  ERROR  FUNCTION  HAVING 

LESS  THAN  1.  3E-7  ERROR  AND  ASSUMED  X  IS  POSITIVE 

ERROF  FUNCTION  IS  SYMMETRIC 

P«0.  327391 1 

A1«0.  234829592 

A2— 0.  284496736 

A3“l.  421413741 

A4— 1.  433132027 

A5«*l.  061405429 

XX-ABS(X) 

T-l.  /(I.  +P*XX> 

ERK»1 .  - ( A1*T+A2*T**2+A3*T**3+A4*T##4+A5*T**3 > *EXP ( -XX**2 ) 

IF  ( X.  OE.  O.  )  ERF“ERF 

IF ( X.  LT.  O.  >  ERF— ERF 

RETURN 

END 
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SUBROUTINE  VTT( JP.  F,  VXTXT,  VXXT.  VTTJ) 

COMMON /ORDER /N, N2 
COMMON/WORNOW/IMODE.  KSMAX.  IOJ 

COMMON/P ARAM/QAMMA < 6. 2).  PHEE(6, 6). H(2,  6). Q<2. 2),  R<2. 2> 
COMMON/P ARAMR /PHEER  (6,6),  DPHEE  (6,6),  GST AR  (6,2),  BSTAR  (2,2) 
DIMENSION  F (6,  6).  VXTXT(6.  6) ,  VXXT(6,  6) ,  VTTJ(2.  2) 

DIMENSION  VHT(6,  2), HVHT(2, 2) , Bi (6, 2) , B2(6, 2). B3(6»  2) 

1,  B4(6,  2),  B5<2,2),  TEMP (6,  6),FL(6.  6) 

DIMENSION  PHEEJ(6,  6),  V(6,  6),  VI  (6,  6),  V2(6,  6),  V3(6,  6),  V4(2.  2) 
IF(JP)  1.2,3 

1  WRITE (6,  11) 

11  FORMAT (2X, 28HNEGAT I VE  J  POWER  IN  VTT  SUB.  ) 

RETURN 

2  CALL  MABCT (H,  2,  N2,  VXTXT,  N2,  H,  2,  HVHT,  2,  6,  6,  6,  2,  6,  2,  2) 

CALL  MATAS(1,  HVHT,  2.  2,  R,  BSTAR. 2, 2) 

DO  4  1  =  1,2 
DO  4  J=l,2 

4  VTTJ(I,  J)=BSTAR(I,  J) 

RETURN 
CONTINUE 

C  V  *  H'  -  GSTAR  *(H*V*H'+R>3 

CALL  MATMUL ( 2,  VXTXT,  N2,  N2,  H,  2,  VHT,  6,  6,  2,  6,  6,  2) 

CALL  MATMUL ( 1 ,  GSTAR,  N2,  2,  BSTAR,  2,  Bl,  6,  2,  2,  2,  6.  2) 

CALL  MATAS (2.  VHT,  N2,  2,  Bl.  B2,  6,  2) 

PHEER  *  C  V  *  H'  -  GSTAR  *(H*V*H'+R)3 
CALL  MATMUL  ( 1,  PHEER,  N2.  N2,  B2,  2,  B3,  6,  6,  6,  2,  6,  2) 

F  =  C  PHEER  *  (  I  -  GSTAR  *  H  )  3**(J-1> 

CALL  CAYLEY ( IMODE.  F, JP-1, FL) 

H  *  F#*(J-1)  *  PHEER  *  C  V  *  H'  -  B  3 

CALL  MATMUL  ( 1 .  FL,  N2.  N2,  B3,  2.  B4,  6,  6.  6,  2,  6,  2) 

CALL  MATMULd,  H,  2.  N2,  B4.  2.  B9.  2,  6,  6,  2.  2,  2) 

IF( IMODE.  EQ.  1)  GO  TO  100 

SUMC  F**(I-1)  *  DPHEE  *  PHEE**(J-I>  1 
DO  5  11  =  1,  N2 
DO  5  12=1, N2 

5  TEMP (II,  I2)=0.  0 
DO  10  1  =  1.  JP 
L1=I— 1 
L2=JP-I 

CALL  CAYLEY( IMODE, F. LI, FL) 

CALL  CAYLEY ( IMODE,  PHEE,  L2,  PHEEJ) 

CALL  MATMULd,  FL.  N2,  N2,  DPHEE,  N2,  VI.  6,  6,  6,  6,  6,  6) 

CALL  MATMULd,  VI,  N2,  N2,  PHEEJ,  N2,  V2,  6,  6,  6,  6,  6,  6) 

CALL  MATAS ( 1 ,  TEMP,  N2,  N2, V2, V, 6, 6 > 

DO  15  11=1, N2 
DO  15  J  J=  1 ,  N2 
TEMP (II,  JJ)=V(II, JJ) 

15  CONTINUE 
10  CONTINUE 

C  SUMC  F**( I— 1 )  #  DPHEE  *  PHEE**(J-I)  3  *  VXXT 

CALL  MATMULd,  V,  N2,  N2,  VXXT,  N2,  V3,  6,  6,  6,  6,  6,  6) 

C  H  *  SUMC  F**(I-1>  *  DPHEE  *  PHEE**(J-I)  3  *  VXXT  *  H' 

CALL  MABCT (H,  2,  N2,  V3,  N2,  H,  2,  V4,  2,  6,  6,  6,  2,  6,  2,  2) 

CALL  MATAS(1,B5,  2,  2,  V4,  VTTJ,  2,  2) 

RETURN 

100  DO  110  1  =  1,2 
DO  110  J=l,  2 
VTTJd,  J)=B5(  I#  J) 

110  CONTINUE 
RETURN 
END 

SUBROUTINE  CAYLEY ( IMODE, F, L, FL) 

C  THI8  SUBROUTINE  PRODUCE  THE  MATRIX  HIGH  POWERED  USING 
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CAYLEY-HAMILTQN'S  THEOREM  TO  REDUCE  I HE  ERROR. 

I  MODE 

(1)  AND  (2)  :  F  IS  DIAOONAL  MATRIX  SOTHAT  IT  HAS  SAME  EIGEN-VALUE. 
THIS  REQUIRE  SPECIAL  OAUS  SUBROUTINE  TO  SOLVE  THE  LINEAR  EQUATIONS. 
(3)  :  F  IS  GENERAL  MATRIX  AND  HAS  THE  DISTINGUISHED 

EIGEN-VLAUE 

F  INPUT  MATRIX  TO  BE  MULTIPLIED  BY  HIGH  POWER 

FL  RESULTANT  MATRIX 

COMMQN/ORDER/N.  N2 

COMPLEX  EV(6), A1 <6, 6).  B1 (6) , ALFA(6> ,  FL1 (6. 6),  CMPLX 
DIMENSION  F(6,  6),  A(12,  12),B(12),  X ( 12) ,  FL(&,  6) ,  FP <6,  6.  6) 

DIMENSION  SF(6,  6) 

IF  (L)  1.2/3 
WRITE (6, 4) 

4  FORMAT ('  NEGATIVE  L  IN  SUB.  CAYLEY') 

RETURN 

DO  9  I-l. N2 

DO  9  J-l.  N2 

FL( I. J)*0. 0 

IF<  I.  EQ.  J)  FL(  I.  J>«1. 

9  CONTINUE 
RETURN 
CONTINUE 

IF(L.  NE.  1 )  GO  TO  7 
DO  6  I-1.N2 
DO  6  J-1.N2 
FL(I.  J)«F(I.  J) 

6  CONTINUE 
RETURN 
CONTINUE 

IFtIMODE.  NE.  3)  GO  TO  190 
CALL  EIGEN(F.  N2.EV) 

USING  GAUSS  ELIMINATION  METHOD,  COMPLEX  MATRIX  CONSISTED  WITH 
EIGENVALUES  IS  PARTITIONED. 

DO  20  I-1.N2 
DO  10  J-1.N2 
10  A1 ( I.  J>— EV( I )**( J-l ) 

20  B 1 ( I ) *EV ( I ) **L 
DO  40  I-1.N2 
DO  30  J-1.N2 
A  ( I • J ) “REAL <  A1 < I ,  J ) > 

Ad.  J^N2)«-AIMA0(A1(I.  J)  ) 

A< I+N2, J)*AIMA0<A1 (I,  J) ) 

Ad  +N2.  J+N2) -REAL (Aid,  J> ) 

30  CONTINUE 

B (I ) -REAL (Bid)) 

B( I+N2)— AIMA0(B1 ( I ) ) 

40  CONTINUE 

CALL  OAUS ( A.  B.  X, N4,  I ERROR) 

GENERATE  THE  COEFFICIENTS  OF  CHARACTERISTIC  FUNCTION 
DO  90  1*1, N2 

ALFA ( I ) -CMPLX ( X ( I > ,  X(I+N2)> 

90  CONTINUE 

CAYLEY-HAMILTON'S  THEOREM 
DO  70  I-1.N2 
DO  70  J-1.N2 
FP (I,  J.  t)»0.  0 
IF(  I.  EQ.  U>  FPd,  J,  1)-1. 

70  CONTINUE 

DO  79  I-1.N2 
DO  79  J-1.N2 
FP(I»J.  2)-F(I,J) 

79  CONTINUE 

NM2-N-2  5? 

IF(NMl)  90,90,99 


95  ucrso  NPS37TT  - - 

DO  85  1  =  1,  N2 
DO  85  J=l,  N2 
FP<  I,  J,  NP  )=0.  0 
DO  85  M=1 , N2 

FP  U .  J. NP  >  =FP U .  J,  NP  ) ♦FP ( I#  M,  NP-1 )*F(H, J) 
85  CONTINUE 
80  CONTINUE 
90  CONTINUE 

DO  100  1*1,  N2 
DO  100  J-1.N2 
FL  1(1,  J)=CMPLX  (O.  O,  0.  0) 

100  CONTINUE 

DO  110  NP=1,  N 
DO  120  1=1, N2 
DO  120  J= 1 , N2 

FL1  ( I,  J)“FL1  <  I,  J >  +ALF A ( NP  )  *FP  ( I >  J,  NP  ) 

120  CONTINUE 
110  CONTINUE 

DO  130  1  =  1.  N2 
DO  130  J=l,  N2 
FLU,  J )  =REAL ( FL  1  ( I ,  J ) ) 

130  CONTINUE 
RETURN 

150  CONTINUE 

DO  152  1=1, N 
DO  152  J=1,N 
152  SF< I.  J)=F(I,  J) 

CALL  EIOEN(SF, N. EV) 

DO  220  1=  1.  N 
DO  210  J=l«  N 
210  A1 < I,  J)=EV( I )**< J— 1 ) 

220  B 1 < I ) =EV ( I ) **L 
DO  240  1=  1.  N 
DO  230  J=1,N 
A( I, J)=REAL<A1< I, J) ) 

A(I,  J+N)=-AIMA0(A1(I,  J) ) 

A( I+N, J ) = A I MAO ( A 1 ( I,  J) ) 

AU+N,  J+N)=REAL(A1<I.  J) ) 

230  CONTINUE 

B ( I ) =REAL ( B 1 ( I >  ) 

B ( I +N ) «A I MAC ( B 1 ( I )  ) 

240  CONTINUE 

CALL  GAUS(A,  B.  X.  N2.  I ERROR) 

DO  250  1=1,  N 

ALFA  < I ) =CMPLX ( X ( I ) ,  X(I+N) ) 

250  CONTINUE 

DO  270  1=  1,  N 

DO  270  J-l.N 

FPU,  J,  1)=0.  0 

IF(  I.  EQ.  J)  FP<  I,  J.  1)  =  1. 

270  CONTINUE 

DO  275  1  =  1,  N 
DO  275  J-l.N 
FPU,  J,  2 )  =SF <  I ,  J) 

275  CONTINUE 
NM2-N-2 

IF(NM2>  290,290,295 
295  DO  280  NP-3, N 
DO  285  I  — 1 ,  N 
DO  285  J-l.N 
FPU,  J.  NP)=0.  0 
DO  285  M-l.N 

FPU.  J.  NP  >=FP  ( I,  J.  NP  )+FPU,  M,  NP  ,i  )*SF(M,  J) 
285  CONTINUE 
280  CONTINUE 


noon  non  nnnn 


290  CONTINUE 

DO  300  I-l.N 
DO  300  J-l.N 
FL1  <  I#  J>— CMPLX(0.  0.0.  0) 

300  CONTINUE 

DO  310  NP-l.N 
DO  320  I-l.N 
DO  320  J-l.N 

FL1  ( I,  J )  — FL  1(1.  J)+ALFA(NP>*FP(I,  J.  NP  > 

320  CONTINUE 
310  CONTINUE 

DO  330  I-l.N 
DO  330  J-l.N 
FL ( I.  J ) —REAL ( FL 1(1.  J)  ) 

FL( I+N.  J+N ) “REAL ( FL 1 ( I . J ) ) 

FL ( I+N. J>-0.  0 
FL(  I.  U+N>=0.  0 
330  CONTINUE 
RETURN 
END 

SUBROUTINE  OAUS(A,  B.  X,  N,  I  ERROR ) 

DIMENSION  A( 12. 12) , B ( 12 > . X ( 12 ) 

C 

C  THIS  SUBROUTINE  IS  IN  'NUMERICAL  ANAYSIS'  BY  L.  W.  JOHNSON  AND  R.  D. 

C  RIESS  , 1977  BY  ADD ISON-WESLEY  PUB.  CO. 

C  SUBROUTINE  OAUS  USES  GAUSS  ELIMINATION  (WITHOUT  PIVOTING)  TO  SOLVE 
C  THE  SYSTEM  AX-B.  THE  CALLING  PROORAM  MUST  SUPPLY  THE  MATRIX  A.  THE 
C  VECTOR  B  AND  AN  INTEGER  N  (WHERE  A  IS  (NXN).  ARRAYS  A  AND  B  ARE 
C  DESTROYED  IN  OAUS.  THE  SOLUTION  IS  RETURNED  IN  X  AND  A  FLAG.  I ERROR. 
C  IS  SET  TO  1  IF  A  IS  NON-SINGULAR  AND  IS  SET  TO  2  IF  A  IS  SINGULAR. 

C  TO  GET  MORE  ACCURATE  SOLUTION.  CALL  SUBROUTINE  RESCOR  AFTER  GAUS. 

C 

NM1-N-1 
DO  9  I-1.NM1 

SEARCH  FOR  NON-ZERO  PIVOT  ELEMENT  AND  INTERCHANGE  ROWS  IF  NECESSARY 
IF  NO  NON-ZERO  PIVOT  ELEMENT  IS  FOUND,  SET  IERROR-2  AND  RETURN 

DO  3  U-I.N 

IF ( A(U>  I).  EQ.  O.  )  GO  TO  3 
DO  2  K-I.N 
TEMP— A ( I,  K) 

Ad.K)-A(U.K) 

2  A(  J.  K )  -TEMP 
TEMP-B ( I ) 

B< I >-B< J) 

B ( J ) -TEMP 
GO  TO  4 

3  CONTINUE 
GO  TO  8 

ELIMINATE  THE  COEFFICIENTS  OF  X ( I )  IN  ROWS  1  +  1.  .  .  .  ,N 

4  IP1-I+1 
DO  9  K-IPl.N 
Q— A ( K,  I) /A(I,  I) 

A(K,  I)-0. O 
B(K)«Q*Bd)+B(K) 

DO  5  J-IP1.N 

9  A  ( K,  J  >  -G*A  ( I ,  J )  +A  ( K»  J ) 

IF(A(N,  N).  EQ.  O.  )  00  TO  S 

BACKSOLVE  THR  EQUIVALENT  TRIANGULAR I ZED  SYSTEM,  SET  IERROR-1. 

AND  RETURN 

X(N)-B(N)/A(N,  N) 
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NPT-N+1 -  - 

DO  7  K-1.NM1 
0*0. 

NMK-N-K 
DO  6  J=l. K 

6  Q«G+A(NMK, NP 1-J ) *X (NP 1-J > 

7  X ( NMK ) * ( B ( NMK ) -Q ) / A  <  NMK,  NMK ) 

I ERROR *1 

RETURN 

8  I ERROR *2 
RETURN 
END 

SUBROUTINE  EIOEN(A, N. EVALUE) 

DIMENSION  A(6, A). RR(6), RI(6), I ANA ( 36 ) >  AT ( 36 > 
COMPLEX  CMPLX, EVALUE(6) 

DO  6  1*1.  N 
DO  7  J*1.N 
K=N* ( I — 1 ) 

7  AT ( J+K ) =A < I »  J) 

6  CONTINUE 

CALL  HSBG(N,  AT.  N) 

CALL  ATEIC(N,  AT,  RR,  RI,  IANA, N) 

DO  5  1  =  1,  N 

EVALUE  ( I  )  =CMPLX  (RRID.RKI)) 

WRITE <6, 500) 

500  FORMAT ( 5X,  'THE  EIGENVALUE  IS') 

WRITE <6, 600) 

600  FORMAT <1 OX.  'REAL  ROOT',  15X,  'I MAG  ROOT') 
WRITE <6, 700)  RR (I).RI(I) 

700  FORMAT  ( 5X.  E15.  6,  14X,  El 5.  6) 

5  CONTINUE 
RETURN 
END 


c 

SUBROUTINE  HSBG 

HSBG 

40 

c 

PURPOSE 

HSBG 

60 

c 

TO  REDUCE  A  REAL  MATRIX  INTO  UPPER  ALMOST  TRIANGULAR  FORM 

HSBG 

70 

c 

USAGE 

HSBG 

90 

c 

CALL  HSBG(N. A,  IA) 

HSBG 

10O 

c 

DESCRIPTION  OF  THE  PARAMETERS 

HSBG 

120 

c 

N  ORDER  OF  THE  MATRIX 

HSBG 

130 

c 

A  THE  INPUT  MATRIX,  N  BY  N 

HSBG 

140 

c 

IA  SIZE  OF  THE  FIRST  DIMENSION  ASSIGNED  TO  THE 

ARRAY 

HSBG 

150 

c 

A  IN  THE  CALLING  PROGRAM  WHEN  THE  MATRIX  IS 

IN 

HSBG 

160 

c 

DOUBLE  SUBSCRIPTED  DATA  STORAGE  MODE.  IA-N 

WHEN 

HSBG 

170 

c 

THE  MATRIX  IS  IN  SSP  VECTOR  STORAGE  MODE. 

HSBG 

180 

c 

HSBG 

190 

c 

REMARKS 

HSBG 

200 

c 

THE  HESSENBERG  FORM  REPLACES  THE  ORIGINAL  MATRIX  IN 

1  THE 

HSBG 

210 

c 

ARRAY  A. 

HSBG 

220 

c 

HSBG 

230 

c 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 

HSBG 

240 

c 

NONE 

HSBG 

250 

c 

HSBG 

260 

c 

METHOD 

HSBG 

270 

c 

SIMILARITY  TRANSFORMATIONS  USING  ELEMENTARY  ELIMINATION 

HSBG 

280 

c 

MATRICES,  WITH  PARTIAL  PIVOTING. 

HSBG 

290 

c 

HSBG 

300 

c 

REFERENCES 

HSBG 

310 

c 

J. H.  WILKINSON  -  THE  ALGEBRAIC  EIGENVALUE  PROBLEM  - 

HSBG 

320 

c 

CLARENDON  PRESS,  OXFORD,  1965. 

HSBG 

330 

c 

HSBG 

340 

c 

HSBG 

350 

c 

HSBG 

360 

SUBROUTINE  HSBG(N, A,  IA) 

HSBG 

370 

DIMENSION  A (36) 

L-N 

55 

HSBG 

400 

n  n  n  n  n  n  non  n  on  non  n  n  n  non 


NIA»L*IA 

LIA-NIA-IA 

L  IS  THE  ROW  INDEX  OF  THE  ELIMINATION 

20  IF(L— 3)  360.40.40 
40  LIA-LIA-IA 
L1«L-1 
L2=L1-1 

SEARCH  FOR  THE  PIVOTAL  ELEMENT  IN  THE  LTH  ROW 

ISUB=LIA+L 
1PIV-IBUB-IA 
PIV*AB8tAt IP IV) ) 

IF(L— 3)  90.90,50 
50  M-IPIV-IA 

DO  80  1*L»  M» IA 
T*ABS(A(I) ) 

IF(T-PIV)  80,80,60 
60  IPIV=I 
PIV=T 

80  CONTINUE 

90  IF(PIV)  100,320,100 
100  IFtPIV— ABSt At ISUB) ) )  180,180,120 

INTERCHANGE  THE  COLUMNS 

120  M-IPIV-L 

DO  140  I*=l »  L 

J«M+I 

T=At J> 

K=LIA+I 
A ( J ) “A ( K ) 

140  A(K)=T 

INTERCHANGE  THE  ROWS 

M=L2-M/IA 

DO  160  I=L1.  NIA.  IA 

T=AtI) 

J=I— M 
At  I )*A< J) 

160  AtJ)«T 

TERMS  OF  THE  ELEMENTARY  TRANSFORMATION 

180  DO  200  I-L, LI A, I A 
200  At  I >*A< I ) /At ISUB  > 

RIGHT  TRANSFORMAT I tm 

J— IA 

DO  240  1-1. L2 

J-U+IA 

LJ-L+U 

DO  220  K-l.Li 

KJ-K+J 

KL-K+LIA 

220  AtKJ)«AtKJ)-AtLJ)»AtKL) 

240  CONTINUE 

LEFT  TRANSFORMATION 

A— IA 

DO  300  I-l.N 


HSBG  410 
HSBG  420 
HSBG  430 
HSBG  440 
HSBG  450 
HSBG  460 
HSBG  470 
HSBG  480 
HSBG  490 
HSBG  500 
HSBG  510 
HSBG  520 
HSBG  530 
HSBG  540 
HSBG  550 
HSBG  560 
HSBG  570 
HSBG  580 
HSBG  590 
HSBG  600 
HSBG  610 
HSBG  620 
HSBG  630 
HSBG  640 
HSBG  650 
HSBG  660 
HSBG  670 
HSBG  680 
HSBG  690 
HSBG  700 
HSBG  710 
HSBG  720 
HSBG  730 
HSBG  740 
HSBG  750 
HSBG  760 
HSBG  770 
HSBG  780 
HSBG  790 
HSBG  800 
HSBG  810 
HSBG  820 
HSBG  830 
HSBG  840 
HSBG  850 
HSBG  860 
HSBG  870 
HSBG  880 
HSBG  890 
HSBG  900 
HSBG  910 
HSBC  920  ' 
HSBG  930 
HSBG  940 
HSBG  950  ' 
HSBQ  960 
HSBG  970 
HSBG  980 
HSBG  990 
HSBG1000 
HSBG1010 
HSBG 1020 
HSBG1030 
HSBG 1040 
HSBG 1050 
HSBG 1060 


K*KVIA - - 

LK=K+L 1 
S-A(LK) 

LJ«=L-IA 
DO  280  J-= 1 , L2 
JK«K+J 
LJ*=LJ-*-IA 

S-S+A<LJ)*A< JK>*1.  ODO 
A  <  LK ) =S 

SET  THE  LOWER  PART  OF  THE  MATRIX  TO  ZERO 

DO  310  I*L> LIA. IA 

A(  I )=0. O 

L=L1 

00  TO  20 

RETURN 

END 


R5BG1070 
HSBC 1 080 
HSB01090 
HSBC 1100 
HSBG1110 
HSBG1120 
HSB01 130 
HSBC 11 40 
HSBQ1150 
HSBC 1 160 
H8BG1170 
HSBG1180 
HSBC 1190 
HSBC 1200 
HSBG1210 
HSBG1220 
HSBG1230 
HSBG1240 
ATE I  10 
.  .  ATE I  20 


SUBROUTINE  ATEIG  ATE I 

ATE  I 

PURPOSE  ATE I 

COMPUTE  THE  EIGENVALUES  OF  A  REAL  ALMOST  TRIANGULAR  MATRIX  ATE  I 

ATE  I 

USAGE  ATE I 

CALL  ATEIG(M, A. RR. RI,  IANA. IA)  ATE  I 


DESCRIPTION  OF  THE  PARAMETERS 
M  ORDER  OF  THE  MATRIX 

A  THE  INPUT  MATRIX.  M  BY  M 


ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 


RR  VECTOR  CONTAINING  THE  REAL  PARTS  OF  THE  EIGENVALUES  ATE I 

ON  RETURN  ATE I 

RI  VECTOR  CONTAINING  THE  IMAGINARY  PARTS  OF  THE  EIGEN-  ATE I 

VALUES  ON  RETURN  ATE I 

IANA  VECTOR  WHOSE  DIMENSION  MUST  BE  GREATER  THAN  OR  EQUAL  ATE I 
TO  M.  CONTAINING  ON  RETURN  INDICATIONS  ABOUT  THE  WAY  ATE I 
THE  EIGENVALUES  APPEARED  (SEE  MATH.  DESCRIPTION)  ATE I 

I A  SIZE  OF  THE  FIRST  DIMENSION  ASSIGNED  TO  THE  ARRAY  A  ATE I 

IN  THE  CALLING  PROGRAM  WHEN  THE  MATRIX  IS  IN  DOUBLE  ATE I 
SUBSCRIPTED  DATA  STORAGE  MODE.  ATE I 

IA=M  WHEN  THE  MATRIX  IS  IN  SSP  VECTOR  STORAGE  MODE.  ATE I 


REMARKS 

THE  ORIGINAL  MATRIX  IS  DESTROYED 

THE  DIMENSION  OF  RR  AND  RI  MUST  BE  GREATER  OR  EQUAL  TO  M 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 

METHOD 

QR  DOUBLE  ITERATION 


ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 
ATE  I 


REFERENCES  ATE I 

J.  O.  F.  FRANCIS  -  THE  QR  TRANSFORMATION - THE  COMPUTER  ATE I 

JOURNAL,  VOL.  4,  NO.  3,  OCTOBER  1961,  VOL.  4,  NO.  4,  JANUARYATEI 
1962.  J.  H.  WILKINSON  -  THE  ALGEBRAIC  EIGENVALUE  PROBLEM  -  ATE I 
CLARENDON  PRESS,  OXFORD,  1965.  ATE I 

ATE  I 


SUBROUTINE  ATEIG (M, A,  RR,  RI,  IANA,  IA) 

DIMENSION  A ( 36 ),RR(6),RI<6),PRR(2),PRI(2), IANA<36> 
INTEGER  P, PI, Q 


ATE  I 
ATE  I 


o  n  n  n  n  n  non  nnnn  non  non 


E7-1.  Ob-8 

AT  El 

490 

E6-1.  OE-6 

ATE  I 

500 

E10-1.  0E-10 

ATE  I 

510 

DELTA-0. 5 

ATE  I 

520 

MAX IT-30 

ATE  I 

530 

ATEI 

540 

INITIALIZATION 

ATE  I 

550 

ATEI 

560 

N-M 

ATEI 

570 

20 

Nl-N-1 

ATEI 

580 

IN— N1*IA 

ATEI 

590 

NN-IN+N 

ATEI 

600 

IF(N1)  30.1300.30 

ATEI 

610 

20 

NP-N+1 

ATEI 

620 

ATEI 

630 

ITERATION  COUNTER 

ATEI 

640 

ATEI 

650 

IT-0 

ATEI 

660  ' 

ATEI 

670 

ROOTS  OF  THE  2ND  ORDER  MAIN  SUBMATRIX  AT  THE 

PREVIOUS 

ATEI 

680 

ITERATION 

ATEI 

690 

ATEI 

700 

DO  40  1-1.2 

ATEI 

710 

PRR(  I )— 0.  0 

ATEI 

720 

40 

PR  I  ( I  >  — 0.  0 

ATEI 

730 

ATEI 

740 

LAST  TWO  SUB DIAGONAL  ELEMENTS  AT  THE  PREVIOUS 

ITERATION 

ATEI 

750 

ATEI 

760 

PAN-0.  0 

ATEI 

770 

PAN1-0. 0 

ATEI 

780 

ATEI 

790 

ORIGIN  SHIFT 

ATEI 

800 

ATEI 

810 

R— 0.  0 

ATEI 

820 

S— 0.  0 

ATEI 

830 

ATEI 

840 

ROOTS  OF  THE  LOWER  MAIN  2  BY  2  SUBMATRIX 

ATEI 

850 

ATEI 

860 

N2-N1-1 

ATEI 

870 

IN1-IN-IA 

ATEI 

880 

NN1-IN1+N 

ATEI 

890 

N1N-IN+N1 

ATEI 

900 

N1N1-IN1+N1 

ATEI 

910 

60 

T— A(N1N1 )— A<NN) 

ATEI 

920 

U-T*T 

ATEI 

930 

V— 4. 0*A ( N 1 N  >  *A ( NN 1 ) 

ATEI 

940 

IF  ( ABS  <  V )  — U*E7 )  100,100.65 

ATEI 

950 

65 

T-U+V 

ATEI 

960 

IF (ABS(T)— AMAX1 (U.  ABS(V>  >*E6>  67. 67. 68 

ATEI 

970 

67 

T— 0.  0 

ATEI 

980 

68 

U«(A(NlNl)+A(NN))/2.  0 

ATEI 

990  , 

V-SQRT(ABS<T))/2.  0 

ATE I 1000 

IF< T )  140.  70.  70 

ATEI 1010 

70 

IF(U)  80.75.75 

ATEI 1020  . 

75 

RR  <N1 >— U+V 

ATEI 1030 

RR ( N ) — U— V 

ATEI 1040 

GO  TO  130 

ATEI 1050 

80 

RR(N1 )— U— V 

ATE I 1060 

RR(N)-U+V 

ATEI 1070 

GO  TO  130 

ATEI 1080 

100 

IF(T) 120, 110, 110 

ATEI 1090 

110 

RR<N1 )— A(N1N1 ) 

ATEI 1100 

RR<N>-A<NN> 

ATEI 11 10 

00  TO  130 

ATEI 1120 

120 

RR(N1 )— A(WO 

58 

ATEI 1130 

RR<N)-A<N1N1) 

ATEI 1140 

non  o  o  o  o  r>  o  r>  ooo 


RI(N1)«0.  0 

ATE Ill 60 

00  TO  160 

ATEI1170 

140 

RR(N1 >*U 

ATEI1180 

RR(N)=U 

ATEI1190 

RI <N1 )*V 

ATE I 1200 

RI(N>=-V 

ATEI1210 

160 

IF(N2) 1280.  1280.  180 

ATEI1220 
ATE I 1230 

TESTS  OF  CONVERGENCE 

ATE I 1240 
ATE I 1250 

180 

N1N2-N1N1-IA 

ATE I 1260 

RM0D=RR(N1 >*RR(N1)+RI(N1)*RI(N1> 

ATE I 1270 

EPS*E 1 0*SQRT ( RMOD ) 

ATE I 1280 

IF( ABS< A(N1N2) )-EPS> 1280.  1280.  240 

ATEI1290 

240 

XF(ABS<A(NN1 ) )-E10*ABS(A(NN) ) )  1300. 1300.250 

ATE I 1300 

250 

IF(ABS(PAN1-A(N1N2) )-AB8(A(NlN2) >*E6)  1240. 1240.  260 

ATEI1310 

260 

IF<ABS(PAN-A<NN1) >-ABS< A(NNl > >*E6 > 1240.  1240. 300 

ATE I 1320 

300 

IF< IT-MAX IT)  320,1240.1240 

ATEI1330 

ATEI1340 

COMPUTE  THE  SHIFT 

ATE I 1350 
ATE I 1360 

320 

J*»1 

ATE I 1370 

DO  360  1=  1.2 

ATE I 1380 

K=NP— I 

ATE I 1390 

IF(ABS(RR<K)-PRR< I ) >+ABS(RI <K>-PRI < I >  >-DELTA*(ABS(RR(K>  > 

ATE I 1400 

1 

+ABS(RI (K) ) > )  340, 360, 360 

ATE I 1410 

340 

U=J+I 

ATEI1420 

360 

CONTINUE 

ATEI1430 

00  TO  (440, 460, 460.  480),  J 

ATE I 1440 

440 

R=0.  0 

ATE I 1450 

S=0.  0 

ATE I 1460 

00  TO  500 

ATEI1470 

460 

J=N+2-J 

ATE I 1480 

R=RR (J) *RR (J) 

ATEI1490 

S=RR(J)+RR(J) 

ATE I 1500 

00  TO  500 

ATE 11510 

480 

R“RR (N)*RR(N1 )— RI (N)*RI (N1 ) 

ATEI1520 

S=RR ( N ) +RR  <  N1 ) 

ATE I 1530 
ATE I 1 540 

SAVE  THE  LAST  TWO  SUBDIAOONAL  TERMS  AND  THE  ROOTS  OF  THE 

ATE I 1550 

SUBMATRIX  BEFORE  ITERATION 

ATE I 1560 
ATE I 1570 

500 

PAN-A(NNl) 

ATE I 1580 

PAN1=A(N1N2> 

ATE I 1590 

DO  520  1=1.2 

ATE I 1600 

K-NP-I 

ATEI1610 

PRR( I )=RR(K) 

ATE I 1620 

520 

PRI(I)«RI(K) 

ATE I 1630 
ATEI1640 

SEARCH  FOR  A  PARTITION  OF  THE  MATRIX,  DEFINED  BY  P  AND  Q 

ATEI1650 
ATE I 1660 

P-N2 

ATEI1670 

IPI-N1N2 

IF  (N-3)  600,600.525 

ATE I 1680 

525 

IPI-N1N2 

DO  580  U-2.N2 

ATEI1690 

IPI«IPI-IA-1 

ATEI1700 

IF(AB8(A< IPI > )-EPS>  600,600,530 

ATEI1710 

530 

IPIP-IPI+IA 

ATE I 1720 

IPIP2-IPIP+IA 

ATEI1730 

D«A( IPIP )*< A( IPIP )-S)+A( IPIP2)*A< IPIP+1 )+R 

ATEI1740 

IF (D) 540,  560.  540 

ATEI1750 

MO  1F(AB8( A< IPX )*A( IPIP+1 ) )*< AB8( A( IPIP )+A( IPIP2+1 >-S)+ABS< A< IPIP2+2) ATE I 1760 
1  >)  -ABS<D)*£PS>  620.620.560  ATEI1770 

560  P-Nl-U  59  ATE X 1 780 


r 

• 

580 

CONTINUE 

ATE I 1790 

600 

Q-P 

ATE I 1800 

00  TO  680 

ATE 11810 

620 

P1=P-1 

ATE I 1820 

a-pi 

IF  (Pl-l)  680,680,650 

650 

DO  660  1*2, PI 

IPI-IPI-IA-1 

ATE I 1850 

IF ( ABS(A( IPI ) ) -EPS > 680,  680,  660 

ATEI1860 

660 

Q=G-1 

ATE I 1870 

C 

ATE I 1880 

C 

OR  DOUBLE  ITERATION 

ATE I 1890 

C 

ATE I 1900 

680 

II=(P-1 )*IA+P 

ATEI1910 

i 

DO  1220  I=P, N1 

ATE I 1920 

I I 1=1 I-IA 

ATE I 1930 

IIP-I1+IA 

ATE I 1940  , 

IF< I-P  >720,  700,  720 

ATE I 1950 

1  700 

IPI-II+1 

ATE I 1960 

IPIP*IIP+1 

ATE I 1970 

C 

ATE I 1980 

C 

INITIALIZATION  OF  THE 

TRANSFORMATION 

ATE I 1990 

C 

ATE 12000 

Q1-A(II)*<A<II)-S>+A<IIP)*A<IPI>+R 

ATE I 20 10 

Q2*A( IPI)*(A( IPIP)+A< II )- 

-s> 

ATE I 2020 

G3*A( IPI )*A< IPIP+1 ) 

ATE I 2030 

A(  IPI  +  1 )=0. 0 

ATE I 2040 

00  TO  780 

ATE 12050 

720 

Q1=A< III) 

ATE I 2060 

G2*A<  1 1 1  +  1 ) 

ATE I 2070 

IF ( I— N2)740» 740,  760 

ATE I 2080 

740 

Q3*A (III +2 ) 

ATE I 2090 

00  TO  780 

ATEI2100 

760 

03*0.  0 

ATE I 21 10 

780 

CAP*SQRT <01*01 +02*02+03*03 ) 

ATE 12120 

IF (CAP >800.  860.  800 

ATEI2130 

800 

IF (01) 820,  840.  840 

ATE 121 40 

820 

CAP— CAP 

ATEI2150 

840 

T-01+CAP 

ATE 121 60 

PSI1-02/T 

ATEI2170 

PSI2*C3/T 

ATE 121 80 

ALPHA*2.  0/ < 1 .  O+PSI 1*PSI 1+PSI2*PSI2> 

ATEI2190 

GO  TO  880 

ATE I 2200 

860 

ALPHA-2.  0 

ATEI2210 

PSI1-0. 0 

ATE I 2220 

PSI2-0. 0 

ATE I 2230 

880 

IF ( 1-0)900,  960.  900 

ATE I 2240 

900 

IFU-P1920.  940.  920 

ATE I 2250 

920 

A<  III)— CAP 

ATE I 2260 

00  TO  960 

ATE 12270 

940 

A  ( 1 1 1 )  —A  (III) 

ATE I 2280  . 

C 

ROW  OPERATION 

ATE I 2300 

c 

ATEI2310 

960 

IJ-II 

ATE I 2320  . 

DO  1040  J«I,N 

ATE I 2330 

T-P8I1*A(IJ+1) 

ATE I 2340 

IF ( I-Nl )980,  1000,  1000 

ATE 12350 

980 

IP2J-IJ+2 

ATE I 2360 

T-T+PSI2*A(IP2J) 

ATE I 2370 

1000 

ETA-ALPHA* ( T+A  < I J) > 

ATE I 2380 

A<  I J)-A< I J)— ETA 

ATE I 2390 

A< I J+l )— A< I J+l >-PSIl*ETA 

ATE I 2400 

IF ( I— N1 ) 1020,  1040,  1040 

ATEI2410 

1020 

A  < IP2J ) -A ( IP2J ) -PS I2*ETA 

ATE I 2420 

1040 

IU-IJ+IA 

60 

ATE I 2430 

,  c 

ATE I 2440 

IB 

•"s' 

non  on  n 


COLUHfrOPCR«TTDW - 

-  -  ATE12450 

ATE I 2460 

IF ( I-Nl ) 1080. 1060.  1060 

ATEI2470 

1060 

K*N 

ATE I 2480 

00  TO  1100 

ATE I 2490 

1080 

K*I+2 

ATE 12500 

1100 

IP*=IIP-I 

ATEI2510 

DO  1180  J=Q, K 

ATE I 2520 

JIP=IP+J 

ATE 12530 

JI=UIP-IA 

ATE I 2540 

T=PSI 1 #A( JIP  > 

ATE I 2550 

IF  < I— N1 ) 1 120.  1140,  1140 

ATE I 2560 

1120 

UIP2<=JIP  +  I  A 

ATE 12570 

T*T+PSI2*A(JIP2> 

ATE I 2580 

1140 

ETA* ALPHA* ( T+A ( J I ) ) 

ATE I 2590 

A(JI)*A(JI>-ETA 

ATE I 2600 

A<UIP)=A(JIP)-ETA*PSI1 

ATE 126 10 

IF ( I— N1 ) 1160,  1180,  1180 

ATE I 2620 

1160 

A(JIP2)»A<JIP2)-ETA*PSI2 

ATE I 2630 

1180 

CONTINUE 

ATE I 2640 

IF< I— N2> 1200,  1220,  1220 

ATE I 2650 

1200 

JI=I 1+3 

ATE I 2660 

UIP=JI+IA 

ATE I 2670 

JIP2=JIP+IA 

ATE I 2680 

ET A-ALPHA*PS I 2* A ( J I P 2 ) 

ATE I 2690 

A( JI )*— ETA 

ATE I 2700 

A ( J I P ) *»-ET A*PS I 1 

ATEI2710 

A(JIP2)*A< JIP2)-ETA*PSI2 

ATE I 2720 

1220 

II-IIP+1 

ATE I 2730 

IT=IT+1 

ATE I 2740 

GO  TO  60 

ATE I 2750 

END  OF  ITERATION 

ATE I 2770 

1240 

IF  <ABS<A<NN1 ) ) — ABS ( A ( N1N2 ) ) ) 

1300. 1280,  1280 

ATE I 2790 

ATE I 2800 

TWO  EIGENVALUES  HAVE  BEEN 

FOUND 

ATEI2810 

ATE I 2820 

1280 

I  ANA  ( N )  *=0 

ATE I 2830 

IANA(N1  )=*2 

ATE I 2840 

N"=N2 

ATE 12850 

IF<N2) 1400,  1400.  20 

ATE I 2860 

ATE I 2870 

ONE  EIGENVALUE  HAS  BEEN  FOUND 

ATE I 2880 

ATE I 2890 

1300 

RR(N)«A(NN) 

ATE I 2900 

RI(N)=0.  0 

ATEI2910 

IANA(N)*1 

ATE I 2920 

IF (N1 ) 1400,  14 JO,  1320 

ATE I 2930 

1320 

N*N1 

ATE I 2940 

GO  TO  20 

ATE 12950 

1400 

RETURN 

ATE I 2960 

END 

ATE I 2970 

61 

m 
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SUBROUTINE  MATMULIIMOT,  A,  N,  M»  B,  L,  C,  NA.  HA,  ' 
DIMENSION  AINA,  MA>,  B(NB,  MB),  CINC,  MC  ) 


A  ,  B  .  C  ARE  GENERAL  MATRIX 

IF  A  X  B  =C.  THEN  IMOT  IS  1 

IF  A  X  B'=C,  THEN  IMOT  IS  2 

DO  1  1*1.  N 

DO  1  J*1.L 

C(I,  J)=0.  0 

DO  1  K»l,  M 

GO  TO  <2, 3), IMOT 

81=BIK, J) 

GO  TO  1 
B1*B( J. K) 

Cl I, J)«CI I.  J)+AU.  K)*B1 

RETURN 

END 

SUBROUTINE  MATAS! I AOS,  A.  N,  M,  B,  C,  NA,  MA> 

DIMENSION  AINA, MA) , B (NA, MA), CINA,  MA) 

IF  A  +  B  =  C,  THEN  IAOS  IS  1 
IF  A  -  B  =  C,  THEN  IAOS  IS  2 
IF!  IAOS.  NE.  1 )  GO  TO  10 
DO  1  1*1, N 
DO  1  J*1,M 

1  C(I,  J)=A(I,  J)+B(I.  J) 

RETURN 

10  DO  2  1*1,  N 
DO  2  J=1,M 

2  C(I.  J>=AII. J)-B ( I,  J) 

RETURN 

END 

SUBROUTINE  MATVECIA,  N. M. B,  C,  NA,  MA) 

DIMENSION  AINA. MA). BIMA), CINA) 

DO  1  1*1,  N 
CII)=0.  0 
DO  1  J-l.M 

C ( I ) *C 1 1 ) +A ( I ,  J)*B(J) 

RETURN 

END 

SUBROUTINE  VECASI IAOS.  A,  B,  C,  N) 

DIMENSION  AIN).  BIN),  CIN) 

A  .  B  ,  C  ARE  VECTORS 
IF  A  +  B  =  C,  THEN  IAOS  IS  1 
IF  A  -  B  *  C,  THEN  IAOS  IS  2 
IF I IAOS.  NE.  1 )G0  TO  10 
DO  1  1*1,  N 

1  C 1 1 >*At I )+BII ) 

RETURN 

10  DO  2  1*1.  N 

2  C 1 1 ) *A 1 1 ) -B 1 1 ) 

RETURN 

END 

SUBROUTINE  MABCTIA,  N,  M,  B,  L,  C,  LL.  D,  NA,  MA, NB, MB, NC, MC, ND. MD) 
DIMENSION  AINA, MA), BINB. MB),  CINC,  MC),  DIND, MD),  ABI6, 6) 

DO  10  1*1. N 
DO  10  U*1,L 
ABII,  J)*0.  O 
DO  10  H-l.M 

ABII,  J)*ABI  I,  J)+AI  I,  K)#B  IK.  J) 

10  CONTINUE 
DO  20  I-l.N 
DO  20  J-l.LL 
DII,  J)*0.  O 
DO  20  K-l.L 

DU.  J)-D!I,  U)-*-ABU,  K)*CI  J.  K) 

20  CONTINUE 
RETURN 
END 


c  ###*#*#**##*#**#**#**#*****#***#*####****###*******###•*##**#* 

C  THIS  PROGRAMMING  IS  CALLED  PHASET.  ITS  MAIN  PURPOSE  IS  TO 

C  ESTIMATE  THE  UNKNOWN  PHASE  USING  THE  PHASE  LOCKED  LOOP  HAVING 

C  A  VERY  NARROW  BANDWIDTH. 

C  PROGRAMMER 

C  CHANG JUNE  YOON 

C  TEXAS  A  *  M  UNIVERSITY 

C  START  JUNE.  1978 

C  #***#*****#***#***##**##*#***#*#******##*###*#**#*#*****•»**** 

COMMON/SAMPLE/NSPB. TB 
COMMON/PHASE/PHEES, PHEEO 
COMMON/GDB/ENODB. SJRDB 

DIMENSION  HMO<2. 2),  HM1<2,  2),  VESTO <4,  4).  VEST1 (4.  4). XESTO(4> 

1. XEST1 (4) . VAR I NO (2,  2 ) >  VAR INI (2.  2) .  VO (2) .  VI (2) 

DIMENSION  GAIN0(4,  2).  GAINK4.  2) 

REAL  MEAN 

LOG I CAL* 1  STRNG  <  8  > 

INTEGER*4  JTIME 

CALL  ASSIGN ( 5.  'SY:  PHASET.  DAT',  13,  'RDO ',  'NC ',  1> 

CALL  INPUT 

READ ( 5,  1 >  NOCASE,  NPRNT 

1  FORMAT <21 5) 

DO  2000  NCASE-1,  NOCASE 
READ <5, 2)  NOSYM, ENODB 

2  FORMAT (15.  E15.  6) 

KSMA  X “NOSYM*NSP  B 

CALL  INITCXJI,  XJQ.  XESTO.  XEST1,  VESTO,  VEST1,  XPI, XPQ, VCO 
1, ERROR, ERRORS.  MEAN. VARANS) 

CALL  OTIM< JTIME) 

CALL  TIMASC< JTIME.  STRNG) 

WRITE<6, 7272)  <STRNC< 1 1 ) .  1 1  =  1,  8) 

7272  FORMAT (IX.  'START  TIME  IS  ',SA1) 

WRITE<6,  50) 

SO  FORMAT <6X. 2HIB.  5X.  5HERR0R.  14X,  6HERRATE.  11X.  6HERR0RS.  12X 
1. 6HERRATS, 12X, 16HPHEE0  IN  DEGREES, 5X. 17HMEAN  AND  VARIANCE) 

DO  JOOO  KS=»1.KSMAX 
CALL  SIONALtKS. BB. SI,  SO) 

CALL  RFKKS.  XJI,  XJQ.  YI.  YQ) 

CALL  DATA<SI.SQ.  YI.  YQ,  ZI,  ZQ> 

CALL  VCOUT <KS,  ZI.  ZQ,  XPI,  XPQ,  VCO,  MEAN,  VARANS) 

CALL  REFOEN(KS»  O.  FTRO,  GTRO,  HMO) 

CALL  REFOEN<KS,  1,  FTR1,  GTR1, HM1 ) 

CALL  KALMAN <KS.  ZI,  ZQ.  HMO, VESTO,  XESTO,  GAINO, VARINO,  DETO, VO) 

CALL  KALMAN (KS,  ZI. ZQ,  HM1, VEST1, XEST1,  GAIN1, VARIN1, DET1, VI ) 

CALL  COSTIKS,  VO.  VARINO,  DETO,  SUMO) 

CALL  COST (KS, VI, VAR INI, DET1, SUM1 ) 

CALL  STAND (K8, ZI.  ZQ,  SUMS,  FTRO,  GTRO,  FTR1, 0TR1 
1,  AFSKO,  AF8K1,  BFSKO,  BFSK1 ,  SFSKO,  SFSK1 ) 

IB-1+IFIX< (KS-. 5) /NSPB ) 

IF(MOD(KS, NSPB ) .  NE.  0)  GO  TO  1000 

CALL  DDC0M(K8, SUMO, SUM1,  XE8T0,  XE8T1, BB, ERROR, ERRATE) 

CALL  STDC0M(K8< SUMS,  SF8K0.  8F8K1,  BB,  ERRORS,  ERRATS) 

PHEE0D-360.  *PHEE0/<2.  *4.  #ATAN( 1.  ) ) 

IF(MOD(IB, NPRNT).  EQ.  0)  WRITE<6,  100)  IB. ERROR, ERRATE, ERRORS,  ERRATS 
1, PHEEOD, MEAN,  VARANS 
100  FORMAT <2X,  15,  5E18.  6.  2E13.  6) 

1000  CONTINUE 

CALL  0TIM< JTIME ) 

CALL  TIMASC (JTIME,  STRNG) 

WRITE (6, 7273)  <STRNQ< II ). II-l.  8) 

7273  FORMAT (IX.  'TIME  IS  ',8A1) 

REWIND  6 

2000  CONTINUE  64 

STOP 


1 

fc.NI) 

C 

BLOCK  DATA 

COMMON/SEED/ IXS, JXS, IXJ1, JXJ1, IXJ2.  JXJ2. IXN1,  JXN1,  IXN2.  JXN2 
COMMON/SAMPLE/NSPB,  TB 
COMMON/OPTION/NOS 
COMMON/DELAY/DELPHI. DELMEG 
COMMON/SIOMA/SICMAJ,  SI OMAN 
COMMON/PHASE/PHEES,  PHEEO 
COMMON /COLOR D/PH I DJ, PHIOJ,  OAMDJ, GAMOJ 
COMMON/ODB/ENODB.  SJRDB 
COMMON/PLLFLT/BNP. EBP. DELF 
COMMON/FREGJ/FJ 

COMMON/PHASIN/HO,  P.  Z,  HI,  KO,  PHA8P.  PHASO 
REAL  KI.KQ 

COMMON/TRACK/OAMMA(4.  4).  PHEE(4, 4) 

INTEGER *2  1X1 (2> . JX1 (2) .  1X2(2), JX2<2>,  1X3(2) , JX3(2) , 1X4(2) . JX4(2> 

1, 1X5(2), JX5(2) 

INTEGER*4  IXS,  JXS,  IXJ1, JXJl.  IXJ2. JXJ2,  IXN1, JXN1,  IXN2,  JXN2 
EQUIVALENCE  (IXS, 1X1), (JXS, JXl), (IXJ1, 1X2), ( JXJ1 , JX2) , (IXJ2, 1X3) 

1, ( JXJ2, JX3) , ( IXN1 , 1X4). ( JXN1. JX4) . (IXN2, 1X5), ( JXN2, JX5) 

DATA  I XI,  JXl/" 136303,  "053354,  "041256,  "141560/ 

DATA  1X2, JX2, 1X3, JX3/"176303,  "037702,  "141236,  "056407. 

1  "125537,  "103453,  "055052,  "032461/ 

DATA  1X4,  JX4,  1X5,  JX5/"034313, "103400,  "021165, "104262, 

1  "072063, "122076. "016415, "041540/ 

END 
C 

SUBROUTINE  INPUT 
COMMON /SAMPLE /NSP B ,  TB 
COMMON/OPTION/NOS 
COMMON/DELAY/DELPHI.  DELMEG 
COMMON/PHASE/PHEES.  PHEEO 
COMMON/FREGJ/FJ 
COMMON/QDB /ENQDB ,  SJRDB 
COMMON/PLLFLT/BNP,  ESP.  DELF 
COMMON/PHASIN/HO,  P,  2,  KI,  KQ,  PHASP,  PHASG 
READ (5,  1)  NSPB.TB 
READ (5,1)  NOS. DELPHI 
READ (5,2)  ENODB . SJRDB 
PHEES=0. 

C  INITIALIZE  PHEES  AS  PHEEO 

PHEE0<=0. 

READ (5, 2)  FJ.  HO 
READ (5, 3)  BNP 

1  FORMAT (15.  El 5.  6) 

2  FORMAT (2E15.  6) 

3  FORMAT  (E15.  6) 

RETURN 
END 

C 

SUBROUTINE  INIT(XJI.  XJQ,  XESTO.  XEST1,  VESTO,  VEST1,  XPI,  XPQ.  VCO 
1, ERROR, ERRORS,  MEAN,  VARANS) 

COMMON/SAMPLE/NSPB,  TB 
COMMON /OP T I ON/NOS 
COMMON/DELAY/DELPHI,  DELMEG 
C0MM0N/8IGMA/S IOMAJ,  SIGMAN 
COMMON/ODB/ENODB,  SJRDB 
COMMON/PLLFLT/BNP,  ESP,  DELF 
COMMON/FREQJ/FJ 

COMMON /COLORD/PH I DJ,  PHIOJ,  OAMDJ,  GAMOJ 
COMMON/PHASE/PHEES,  PHEEO 
COMMON/PHASIN/HO, P, Z, KI,  KQ,  PHASP, PHASO 
REAL  KI.KQ,  MEAN 

COMMON/TRACK/OAMMAU,  4),  PHEE(4,  4) 

DIMENSION  XESTO ( 4 ) , XEST 1(4), VESTO (4,4). VE8T 1(4,4) 


65 


non  nnnn  o  nn 


r - - - 

PI-4.  *  AT  AN  ( 1 .  ) 

DELMEO-DELPH I *2 .  *PI/TB 
IF(NOS.  EQ.  1 )  00  TO  10 
SUMF-0. 

DO  15  K-l.NSPB 

15  SUMF— SUMF+(SIN< (K-.  5)*TB*DELME0/NSPB) )**2 
SI OMAN— SORT  <SUMF)*10.  **( -EN0DB/20.  > 

00  TO  20 

10  CONSTP— SORT (NSPB/2.  )*ABS(SIN(DELPHI ) ) 

SI0MAN-C0NSTP*10.  **<-EN0DB/20.  ) 

20  SIGMAJ— 10.  **(— SJRDB/20.  >/SQRT(2.  > 

GENERATE  THE  COLOURED  NOISE  PARAMETERS  AND  ITS  BANDWIDTH 
"  THE  RHO-FILTER  AND  ITS  BANDWIDTH 
T-TB/NSPB 
POLE J— 2.  *PI*FJ 
PHIDJ-EXP  <POLEJ*T ) 

CAM— (PHI D J- 1 .  ) /POLEJ 
OAINK-1.  /SORT ( GAM#*2/ ( 1 .  -PHIDJ**2> ) 

0 AMD J— 0 A I NK*GAM 
PHIOJ-O. 

OAMOJ-O. 

BNJ— POLEJ/4. 

BNR-BNP 
P0LER--4.  *BNR 
PH I DR-EXP (POLER*T ) 

CAM— (PHIDR— 1 .  ) /POLER 
OAINK-1.  /SQRT(0AM**2/(1.  -PHIDR**2> ) 

GAMDR-OA I NK*QAM 
PHIOR-O. 

OAMOR-O. 

DO  50  1-1.4 
DO  50  J-l. 4 
GAMMA ( I.  J)— O. 

50  PHEE<  I.  J)— O. 

GAMMA  ( 1 .  D-GAMDR 
GAMMA (2.2) -GAMDR 
GAMMA  <  3. 3 ) =0 AMD J 
GAMMA (4, 4  >  — CAMD J 
PHEEtl. 1 ) -PHIDR 
PHEE<2. 2) -PHIDR 
PHEE(3.  3)— PHIDJ 
PHEE<4.  4>— PHIDJ 

GENERATE  THE  PHASE  ESTIMATOR  PARAMETERS 
A— 2.  *PI*ESP/360. 

TANHO-SIN ( A ) /COS  <  A ) 

HO- ( 2.  *P I*DELF ) /TANHO 
KG- (8.  /3.  )*BNP 
Z—  <4.  /3.  )*BNP 
P-KG*Z/HO 
KI-P/Z 

PHASP— EXP  <  P*T ) 

PHASO— ( PHASP— 1 .  )/P 

INITIALIZATION 


XJI-O. 

XJG-O. 

XPQ-O. 

XPI-1.  / <KI*(P-Z ) ) 
VCO-O. 

DO  60  1-1.4 
XE3T0( I )— 0. 

60  XEST1  ( I )— 0. 

DO  65  1-1,4 
DO  65  J-l.  4 


L 


^ --fi8  «i  fn  'n 


AJk 


VfcSHM  1.  Jl«U. 

IFd.EQ.  J)  VESTOd.  J)-l. 

65  CONTINUE 
DO  70  1-1.4 
DO  70  J=1 , 4 

70  VEST1 ( I.  J)=VESTO< I.  J) 

ERROR-O. 

ERRORS— 0. 

PHEE0=0. 

MEAN-0. 

VARANS— O. 

WRITE <6. 99)  ENODB.SJRDB 

99  FORMAT (2X,  6HEN0DB-,  E13.  6.  5X,  6HSJRDB-,  E13.  6.  /) 

WRITE (6,  100)  NOS.  NSPB.  TB.  DELPHI. PHEES 

100  FORMAT <2X.  4HN0S-.  12.  5X.  5HNSPB-.  15.  5X. 3HTB-,  E13.  6.  5X 
1.  7HDELPHI-.  E13.  6.  5X.  6HPHEES-.  E13.  6.  /> 

WRITE (6. 101 )  GAMDJ.  PHIDJ.  BNJ 

101  FORMAT ( 5X.  6HGAMDJ-. E13.  6. 5X.  6HPHIDJ-,  E13.  6. 5X,  4HBNJ-.  El 3.  6) 

WRITE (6,  102)  GAMDR, PHIDR.  BNR 

102  FORMAT <5X.  6HGAMDR-,  E13.  6.  SX.  6HPHIDR-.  E13.  6.  5X.  4HBNR-.  E13.  6) 

WRITE <6,  103)  PHASG.  PHASP, BNP 

103  FORMAT <5X,  6HPHASG-,  E13.  6.  5X.  6HPHASP-.  E13.  6.  SX.  4HBNP-.  E13.  6) 
WRITE<6,  105)  HO.  P.Z.KI.KQ 

105  FORMAT (2X.  18HPARAMETERS  IN  VCO-, />  5X<  SHH(0)>>  E13.  6*  5X.  2HP-.  E13.  6 
1.  5X<  2HZ— .  E13.  6,  5X.  3HKI-.  E13.  6.  SX.  3HKQ-.  E13.  6.  ///> 

REWIND  6 

RETURN 

END 

SUBROUTINE  SIGNAL(KS.  BB.  SI.  SQ) 

COMMON/SEED/ 1 XS.  JXS,  IXJ1.  JXJ1,  IXJ2,  JXJ2.  IXN1.  JXN1.  IXN2.  JXN2 

INTEGER *4  IXS, JXS.  IXJ1.  JXJ1.  IXJ2,  JXJ2.  IXN1,  JXN1.  IXN2.  JXN2 

COMMON/SAMPLE/NSPB. TB 

COMMON/OPTION/NOS 

COMMON/PHASE /PHEES. PHEEO 

COMMON/DELAY/DELPHI.  DELMEG 

IF(M0D(KS-1, NSPB).  NE.  0)  GO  TO  10 

CALL  RANCdXS,  JXS,  QB) 

BB-AINT(GB+.  5) 

10  C-l.  — 2*BB 

TK— (KS-.  5) /NSPB 
TKMOD- ( TK-IFI X ( TK ) )*TB 
A-l. 

GO  TO  (1,2),  NOS 

1  PHEEM-DELPH I *C 
GO  TO  20 

2  PHEEM-DELMEC*C*TKMOD 
20  SI— A*COS ( PHEEM+PHEES > 

SQ-A*S I N ( PHEEM+PHEES ) 

RETURN 

END 

SUBROUTINE  RFKKS,  XJI,  XJG,  YI,  YQ) 

COMMON/SEED/ IXS,  JXS,  IXJ1, JXJ1, IXJ2. JXJ2, IXN1, JXN1,  IXN2, JXN2 
INTEGER *4  IXS.  JXS, IXJ1.  JXJ1.  IXJ2, JXJ2. IXN1, JXN1,  IXN2,  JXN2 
COMMON/COLORD/PH I DJ,  PHIOJ.  GAMDJ,  OAMOJ 
COMMON/S I OMA/S I OMAJ,  SIOMAN 
REAL  NI.NQ 

CALL  MARSAdXJl,  JXJ1.  WI ) 

CALL  MARSAC IXJ2,  JXJ2.  WQ) 

CALL  MARSAdXNl,  JXN1,  NI  > 

CALL  MARSA( IXN2,  JXN2,  NG) 

X J I 1 -PH I D J*X J I +PHI 0 J*X JQ+OAMD J*W I +OAMO J*WO 

X JO 1 —PH I 0 J#X J I +PH I D J*X JG-0 AMO J*W I +OAMDJ*WQ  6 

YI»SIOMAJ*XJI+SIOMAN*NI 


yq-s igmau*  x ju+s 1 gman*nq 

XJI-XJU 

XJQ-XJQ1 

RETURN 

END 

C 

SUBROUTINE  DATA ( SI >  SG,  YI»  YG,  ZI,  ZG) 

COMMON/PHASE/PHEES, PHEEO 

ZI-SI+YI 

ZG-SG+YG 

Z I  «Z  I *COS ( PHEEO  >  +ZG*S IN ( PHEEO ) 

ZQ--Z I*S IN ( PHEEO ) +ZG*COS( PHEEO) 

RETURN 

END 

C 

SUBROUTINE  VCOUT(KS,  ZI, ZQ,  XPI,  XPQ,  VCO.  MEAN/  VARANS) 

COMMON/SAMPLE/NSPB,  TB 

COMMON/PHASE/PHEES, PHEEO 

COMMON/PHASIN/HO.  P,  Z,  KI,  KG,  PHASP.  PHASG 

REAL  KI. KG,  MEAN 

XPQ1«PHASP*XPQ+PHASG*ZQ 

ZQ1-KQ*( (P-Z)*XPQ+ZG) 

XPG-XPG1 

XP I i»PHASP*XP I+PHASG*Z I 
ZI1-KI*((P-Z)*XPI+ZI) 

XPI-XPI1 
VCQP1-ZG1/ZI1 
TCTB /NSPB 

PHEEO—PHEEO+  ( VCOP 1 +VCO )  *T/2. 

VCO=VCOP 1 

C  ESTIMATE  THE  MEAN  AND  VARIANCE  OF  THE  PHASE  ERROR.  RECURSIVELY 

PHEE0D-360.  *PHEEO/ (2.  *4.  *AtAN( 1.  ) ) 

MEAN— ( (KS— 1.  )*MEAN+PHEEOD)/KS 
EVAR- ( P HEEOD— ME AN ) **2 
VARANS- ( <  KS- 1 .  ) *VAR ANS+E VAR ) /KS 
RETURN 
END 
C 

SUBROUTINE  REFGEN(KS.  M,  FTR,  GTR. HM) 

C0MM0N/8AMPLE/NSPB, TB 
COMMON/DELAY/DELPHI,  DELMEG 
COMMON/OPTION/NOS 
DIMENSION  HM(2,  2) 

TK— (KS—.  5) /NSPB 
TKMOD- ( TK- IF I X ( TK ) ) *TB 
AR-1. 

IF (NOS.  NE.  1)  00  TO  1 
IF(M.  EG.  O)  PHEEMR-DELPHI 
IF(M.  EG.  1)  PHEEMR-DELPHI 
GO  TO  2 

1  IF(M.  EG.  0)  PHEEMR«DELMEC#TKMOD 
IF(M.  EG.  1)  PHEEMR— -DELMEQ*TKMOD 

2  FTR«AR*COS( PHEEMR) 

QTR»AR*S IN (PHEEMR) 

HM( 1, 1 > -COS (PHEEMR) 

HM( 1,  2>— SIN(PHEEMR ) 

HM(2.  l)-SIN(PHEEMR) 

HM(2.  2)— COS  (PHEEMR) 

RETURN 

END 

C 

SUBROUTINE  KALMAN (KS,  ZI,  ZQ,  HM.  VEST,  XEST, GAIN, VARINV, DET,  V) 
COMMON/TRACK/GAMMA (4,  4),  PHEEC4,  4) 

COMMON/SIGMA/SIGMAJ,  SIQMAN 

DIMENSION  VEST (4,4),  PVP ( 4,  4 > ,  GTG (4,4),  VPRED ( 4,  4),  VHT(4,  2) 

1,  HVHT(2,  2) ,  VAR (2,  2),  VARINV(2.  2),  GAIN(4,  2),  GH(4,  4),  HM(2,  2) 


mwertSrDRT  VNN(2.  2) 

REAL  IMGH(4, 4) 

DIMENSION  XE8T ( 4 ) > XPRED(4>,  HXPRED(2>, V(2>, GV(4>,  HX(2,  4> 
DO  1  1-1.2 
DO  1  J-1,2 
1  HXd,  J)«HM(I,  J) 

HX ( 1 < 3)-SIGMAJ 
HX ( 1 .  4>-0. 

HX  (2.  3)— 0. 

HX(2> 4>— SIOMAJ 
VNN(1, 1 >— SI0MAN**2 
VNN(1,  2)-0. 

VNN(2,  l>-0. 

VNN<2.  2 ) — S 1 CMAN**2 

C  CALCULATE  THE  STEADY-STATE  KALMAN  GAIN 

CALL  MABCT (PHEE.  4,  4.  VEST.  4.  PHEE,  4.  PVP.  4.  4.  4.  4.  4.  4,  4.  4) 
CALL  MATMUL ( 2.  GAMMA.  4.  4.  GAMMA.  4.  OTG.  4.  4.  4.  4.  4.  4) 

CALL  MATAS< 1.  PVP, 4. 4,  OTG, VPRED, 4.  4) 

CALL  MABCT <HX,  2,  4.  VPRED.  4,  HX.  2,  HVHT,  2,  4,  4,  4.  2,  4.  2.  2> 
CALL  MATAS ( 1 .  HVHT,  2,  2.  VNN.  VAR,  2,  2  > 

DET-VARd,  1>*VAR(2,  2)-VAR(l,  2)*VAR(2,  1) 

VARINVd,  1  >  -VAR  (2,2)  /DET 
VARINVd,  2 )  —VAR ( 1 ,  2) /DET 
VARINV(2> 1 )=-VAR(2» 1)/DET 
VARINV (2,2) -VAR ( 1 . 1)/DET 

CALL  MATMUL < 2.  VPRED.  4.  4,  HX,  2,  VHT ,  4,  4,  2,  4,  4.  2) 

CALL  MATMUL (1.  VHT,  4,  2,  VARINV.  2,  GAIN,  4,  2,  2,  2,  4,  2) 

CALL  MATMUL ( 1.  GAIN,  4.  2.  HX.  4.  OH,  4,  2,  2,  4,  4,  4) 

DO  10  1-1.4 
DO  10  J-1,4 
IMGHd,  J)  — GHd,  J) 

IFd.EQ.  J)  IMGHd,  J>-1.  -GHd,  J) 

10  CONTINUE 

CALL  MATMUL < 1 ,  IMOH,  4.  4,  VPRED, 4. VEST. 4. 4, 4, 4. 4. 4) 

C 

CALL  MATVEC (PHEE. 4. 4.  XEST,  XPRED. 4, 4) 

CALL  MATVEC < HX ,  2.  4,  XPRED,  HXPRED,  2,  4) 

V(l)-ZI-HXPREDd) 

V ( 2 ) -ZG— HXPRED  <  2 ) 

CALL  MATVEC (GAIN, 4,  2.  V,  OV,  4.  2) 

CALL  VEC AS ( 1 , XPRED.  GV,  XEST.  4) 

RETURN 

END 

C 

SUBROUTINE  COST(KS,  V,  VARINV,  DET, SUM) 

COMMON /SAMPLE /NSPB ,  TB 
DIMENSION  V(2)>  VARINV (2.  2) 

IF(M0D(KS-1»  NSPB).  EQ.  0)  SUM-0. 

ARO— ALOO  ( DET )  -  ( V  ( 1 )  **2*VAR  INV  d ,  1 )  +V  ( 2 )  **2*VAR  I NV  ( 2.  2  > 
1+V( 1 )*V(2)*( VARINV( 1,  2)+VARINV(2,  1 ) ) ) 

3UM-SUM+ARG 
RETURN 
END 
C 

SUBROUTINE  STAND (KS,  ZI.  ZQ,  SUM,  FTRO.  GTRO,  FTR1,  GTR1 
1, AFSKO, AFSK1,  BFSKO.  BF8K1,  SFSKO, SFSK1 > 

COMMON /SAMPLE /NSPB.  TB 
COMMON/OPTION/NOS 
GO  TO  (1, 2), NOS 

1  IF(M0D(KS-1,  NSPB).  EQ.  0)  SUM-0. 

SUM— SUM+ZQ 
RETURN 

2  IF(M0D(KS-1, NSPB).  NE.  0)  GO  TO  20 
AFSKO— 0. 

AFSK1-0. 

BFSKO— 0. 


SHE 


Bt-bni-O. 

20  AF8K0»AFBK0+ZI*FTR0+ZQ*©TR0 
AF8K 1  -AFSK  1 +Z I  *FTR  1  +ZQ*OTR  1 
BFSKO-BFSKO+ZI*OTRO-ZQ*FTRO 
*  BF8K 1 -BF8K 1 +Z I *©TR 1 -ZQ*FTR 1 
IF(MOD(KS,  NSPB ) .  NE.  0)  RETURN 
SFSK0«AFSK0**2+BFSK0**2 
SFSK 1 -AFSK 1 **2+BFSK 1 **2 
RETURN 
END 

SUBROUTINE  DDCOMWS,  SUMO,  SUM1.  XESTO,  XEST1,  BB,  ERROR,  ERRATE) 
COMMON/BAMPLE/NSPB, TB 
DIMENSION  XESTO ( 4 ) , XEST 1(4) 

IF (SUMO.  OT.  SUM1)  GO  TO  10 
BBHAT- 1. 

DO  1  1-1,4 

1  XESTO< I )-XESTl ( I ) 

00  TO  20 

10  BBHAT-O. 

DO  2  1-1,4 

2  XEST1 < I )— XE8T0( I > 

20  IF<BB.  EQ.  BBHAT)  ERR-O. 

IF(BB.  NE.  BBHAT)  ERR-1. 

ERROR-ERROR+ERR 

IB— 1+IFIX( (KS-.  5) /NSPB) 

ERR ATE-ERROR /IB 

RETURN 

END 

SUBROUTINE  STDCOM(KS.  SUM,  SFSKO,  8F8K1.  BB. ERROR, ERRATE) 
COMMON/SAMPLE/NSPB.  TB 
COMMON/OPTION/NOS 
00  TO  (1.2).  NOS 

1  IF(SUM.  OE.  O.  )  BBHAT-O. 

IF  (SUM.  LT.  0.  )  BBHAT- 1. 

00  TO  10 

2  IF ( SFSKO.  OT.  SFSK 1)  BBHAT-O. 

IF(SFSK1.  OT.  SFSKO)  BBHAT-1. 

10  IF ( BB.  EQ.  BBHAT)  ERR-O. 

IF ( BB.  NE.  BBHAT)  ERR-1. 

IB— 1+IFIX( (KS-.  S ) /NSPB > 

ERROR-ERROR+ERR 
ERR ATE-ERROR /IB 
RETURN 
END 

SUBROUTINE  MARSA( IXA.  JXA.  V) 

INTE0ER*4  IXA,  JXA 
CALL  RANCdXA,  JXA,  XI) 

CALL  RANCdXA.  JXA.  X2> 

X1»(X1-.  9)  *2. 

X2— (X2-.  5) *2. 

9  M-X1**2+X2**2 

IF(W.  LE.  1.  )  00  TO  10 
CALL  RANCdXA.  UXA.  XI) 

CALL  RANCdXA,  UXA.  X2> 

X1-(X1-.  9) *2. 

X2-(  X2-.  9)  *2. 

00  TO  9 

10  XX-Xl*SGRT(-2.  *ALOO(W)/W) 

V— X2*XX/X1 

RETURN 

END 
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